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ABSTRACT 

(Distribution  Limitation  Statement  No.  2) 


An  integrated  computer  program  (STRAB-6)  for  the  complete  analysis  of  a  blunt, 
conical  reentry  vehicle  as  to  its  trajectory  and  aerothermal  environment  upon 
atmospheric  reentry  is  presented.  The  trajectory  portion  of  the  program  calcu¬ 
lates  the  vehicle  notions  in  six-degrees-of-freedom,  while  the  thermal  portion 
calculates  the  aerodynamic  heating  and  ablation  of  the  vehicle.  The  earth 
model  is  an  oblate,  rotating  spheroid  whose  parameters  and  gravitational 
potential  are  described  herein.  The  equations  of  motion,  the  method  of  integra¬ 
tion,  thermal  model  and  input  forms  for  the  program  are  discussed.  The  thermal 
model  examined  in  this  report  assumes  thick  skin  solution  where  the  skin  is  of 
any  composite  structure  made  of  discrete  lavers  of  material  whose  properties 
may  vary  from  layer  to  layer.  Also,  the  heatshield  is  comprised  of  one,  two, 
or  more  different  ablative  materials.  STRAB-6  computes  nose-blunting  and 
includes  this  effect  in  the  aerodynamics.  STRAB-6  is  written  in  FORTRAN  IV 
for  the  CDC  6600  computer. 
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SECTION*  I 
INTRODUCTION 


Tne  effects  of  aerodynamic  heating  ar.d  ablation  on  rhe  motions,  aerodynamic 
forces,  and  trajectory  deviations  of  a  reentry  vehicle  have  been  the  subject  of 
cany  analyses  in  recent  years-  Previous  trajectory  and  heating  programs  have 
been  separate,  requiring  that  each  be  run  separate1 y  with  data  from  one  used  as 
input  to  the  other.  Such  a  method  sometimes  requires  numerous  runs  before  the 
output  data  converge  to  reasonable  answers,  and  is  both  expensive  in  computer 
costs  and  time  consuming.  This  report  presents  a  computer  program  (STKAB-6) 
for  the  cccplete  analysis  of  a  reentry  vehicle  as  to  its  trajectory  and  aero- 
thermal  environment.  STRAB-6  provides  a  six-degree-of-freedon  trajectory 
analysis  as  well  as  aerodynamic  heating  and  ablation  analysis  for  conical 
vehicles  reentering  the  atmosphere  at  hypersonic  speeds. 

STRAB-6  uses  seme  of  the  features  from  other  6-D  programs  (Refs-  1,  2,  3, 

4,  5,  6).  The  trajectory  portion  of  the  program  computes  the  six-degree-cf- 
freedom  rigid  body  motions  for  an  aerodynamical ly  symmetric  reentry  vehicle 
encountering  a  model  atmosphere.  Slender-body  aerodynamic  theory  was  used  and 
the  required  aerodynamic  coefficients  are  provided  by  subroutine  input-*  (as  a 
function  vehicle  shape,  angle-of-attack,  and  either  Mach  number  or  altitude). 
The  atmospheric  model  that  is  used  is  based  on  the  1962  ARDC  atmosphere,  and 
the  geocentric  model  that  is  used  is  a  rotating  oblate  spheroid. 

Tne  aerodynamic  heating  and  ablation  portion  of  the  program  (Kef.  7) 
(HFATAB),  in  its  original  form,  computed  point  mass  heating  for  a  sharp-ncsec 
conical  vehicle.  The  basic  heating  equations  of  this  original  program  (KEA7AE) 
are  described  in  reference  7.  There  have  been  numerous  improvements  to  the 
original  program  to  include:  heating  on  a  blunt-nosed  vehicle,  the  case  where 
a  vehicle  he3tshield  is  comprised  of  two  or  more  different  ablators,  and 
boundary  layer  edge  conditions. 


The  equations  of  motion  and  of  heat  transfer  are  general  in  that 
to  most  axially  symmetric  vehicles.  Even  though  this  program  conside 
body  in  general  motion  under  combined  angles  of  attack  and  sideslip, 
vehicle  is  assumed  to  ablate  symmetrically.  The  ablation  weight  lc-s*- 
bluntness,  and  vehicle  volume  changes  are  considered  ir.  the  calculati 
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the  viscid  and  ir.?iscid  drag  coefficients  (CD)  and  stability  derivatives,  as 
well  as  the  vehicle  mass  moment  of  inertia  properties.  In  addition,  the  aero¬ 
dynamic  coefficients  are  calculated  for  the  complete  range  of  angle-of-attack 
(0  <  a  <  360°)  and  sideslip  (0  ^  S  360°). 

The  position  vector  is  calculated  in  a  rotating  esrtn-centered  coordinate 
system,  while  the  angular  positions  are  calculated  in  a  reentry-centered 
coordinate  system.  Computations  are  carried  out,  using  the  Runge-Kutta  method 
of  numerical  integration  employing  a  fired  step  size-mode. 

This  program  is  specific  in  that  the  equations  have  been  applied  to  a 
blunt-nosed  conical  vehicle  with  the  starting  point  at  reentry  conditions. 

With  some  edifications ,  the  user  could  easily  adapt  the  program  to  analyze  a 
missile  from  launch  by  adding  thrust  subroutines  and  substituting  launch 
conditions  for  reentry  conditions- 
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SECTION  II 

MAIN  PROGRAM  DESCRIPTION 

1.  GENERAL 

The  computer  program,  STRA3-6,  is  a  self-contained  programs  in  the  sense 
that  it  can  generate  the  required  aerodynamic  data,  and  calculate  the  mass  and 
inertial  properties  due  to  ablation  changes  of  a  symmetrical  vehicle  without 
any  additional  preliminary  calculations. 

STRAS-6  calculates  the  resultant  loads,  motions,  and  trajectory  of  a  rigid 
vehicle  in  general  6-D  notion  as  veil  as  boundary-layer  properties  and  time- 
temperature  distribution  in  a  body-  If  a  conablating  trajectory  is  desired, 
a  option  is  provided  to  omit  the  heating  and  ablation  portion  of  the  program. 

At  the  initiation  of  a  reentry  run,  the  vehicle  is  located  in  space  with 
the  desired  reentry  conditions  (i.e.,  velocity,  flight  path  angle,  altitude, 
roll  rate).  Also,  the  user  must  input  the  latitude,  longitude,  and  azimuth  at 
which  he  wishes  to  start  the  trajectory.  Other  trajectory  parameters  the  user 
nay  vish  to  input  are  angle— of-attack  and  angle— of-sideslip;  otherwise,  these 
values  are  zero  for  a  nominal  trajectory. 

The  aerodynamic  drag  coefficients  (CD)  are  obtained  from  Newtonian  siender- 
body  theory  for  inviscid  drag  and  modified  Bias Ins  flat  plate  theory  for 
viscous  crag-  These  coefficients  are  calculated  in  subroutines  CDS?  and  CD?X 
and  are  limited  to  vehicle  shapes  of  spherically  blunted  cones  where  the  ratio 
of  base  radius  to  axial  length  is  approximately  equal  to  the  cone  half-angle. 

A  discussion  of  these  theories  is  in  section  I 1-3. 

The  atmospheric  model  used  In  the  urogram  Is  based  on  the  1962  AEDC  model 
atmosphere,  and  it  expresses  the  density,  ambient  temperature,  pressure,  and 
the  speed  of  sound  as  functions  of  altitude.  The  geocentric  model  -used  is  a 
rotating  oblate  spheroid,  where  the  sea  level  radius,  S_,  is  (figure  1) 
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where 


R£p  =  earth  equa.orial  radius-ft 
=  geocentric  latitude-degree 


P 


E 


=  0,00673852  (Ref  8) 


e£  =  eccentricity  of  earth  =  0.0818133302 


(2) 


Figure  1.  Geocentric  Earth  Model 

The  vehicle  thermal  model  as  presented  in  HEATAB  has  been  modified  to 
extend  the  calculations  to  include  the  total  integrated  ablation  mass  loss  rate 
for  a  conical  reentry  vehicle.  The  body  is  divided  into  a  number  of  segments, 
or  stations  along  the  cone,  and  the  heating  rate  at  each  station  is  computed. 

In  addition  to  the  original  input  data  required  by  HEATAB,  the  fore  and  aft 
radii,  slant  length,  and  transition  Reynolds  number  of  .ach  body  segment  must 
be  supplied.  These  input  data  are  identified  as:  Data  RS,  Data  RL,  Data  SOL, 
and  Data  RN.  Computations  of  the  skin  temperature  profile  and  ablation  mass 
loss  rate  are  identical  to  the  method  in  HEATAB,  except  that  these  calculations 
are  made  for  each  X-station  at  each  trajectory  point. 
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As  written,  the  STRAB-6  output  will  give,  for  each  time  interval,  an  output 
of  trajectory  conditions,  body  conditions,  and  aerodynamic  parameters.  The 
trajectory  output  conditions  include  time,  altitude,  velocity,  and  flight  path 
angle.  The  body  conditions  for  each  X-station  include  weight  loss  rate  per 
unit  area,  wall  recession,  and  wall  temperature  profile.  The  aerodynamic 
parameters  include  CD,  maximum  aerodynamic  load,  angle-of-attack  and  sideslip, 
restoring  forces  and  moments,  and  boundary  layer  edge  conditions. 

2.  EQUATIONS  OF  MOTION 

The  equations  of  motion  for  an  ascent  or  descent  through  the  atmosphere 
can  be  written  with  various  degrees  of  complexity.  The  equations  differ  in 
form  depending  upon  the  coordinate  system,  the  number  of  degrees  of  freedom, 

the  method  of  adjusting  for  the  earth’s  rotation,  and  the  integration  scheme. 

< 

In  the  evaluation  of  reentry  vehicles  for  height-of-burst  (HOB)  errors,  disper¬ 
sion  analysis,  or  circular  error  probable  (CEP),  the  equations  of  morion  must 
be  written  with  six-degrees-of-freedom  (6-D).  This  program  employs  a  Runge- 
Kutta  method  of  integration  for  a  high  degree  of  accuracy. 

a.  Coordinate  Systems 

The  equations  of  motion  can  be  written  in  vector  notation  using  Newton's 
second  law.  A  summation  of  the  forces  acting  on  the  vehicle  results  in  the 
equation 


Before  this  equation  can  be  solved,  however,  it  must  be  expressed  in  terms  of 
a  convenient  coordinate  system. 

There  are  five  principal  coordinate  systems  or  reference  frames  associ¬ 
ated  with  orbital  mechanics  and  trajectories.  Various  other  reference  frames 
can  be  used  depending  on  the  complexity  of  the  problem  to  be  solved.  For  this 
report,  there  are  three  reference  frames:  inertial  reference  frame,  topocentric 
or  horizontal  reference  frame,  and  earth  reference  frame.  The  following  is  a 
description  of  each  and  their  interrelationships.  All  the  reference  frames 
are  right-handed  orthogonal  coordinate  systems. 
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(1)  Inertial  Reference  Frame 

The  basic  coordinate  for  the  program  is  a  nonrotating  inertial 
reference  frame  fixed  with  its  origin  at  the  geocenter  as  shown  in  figure  2. 

ZR  points  north  along  the  axis  of  rotation  of  the  earth  and  XR  and  YR  lie  in 
the  equatorial  plane.  XR  lies  in  the  prime  meridian  plane  at  the  start  of  a 
trajectory  run  (time,  t  =■  0).  The  earth  rotates  with  an  angular  velocity  w£ 
about  ZR;  hence  the  earth  reference  frame,  E,  differs  from  the  inertial  frame, 

R,  by  a  rotation  w£t. 

The  vehicle  center  of  gravity  (CG)  is  located  in  the  inertial  R 
frame  by  the  longitude,  8,  latitude,  Xc,  and  vehicle  altitude,  ALT.  The 
distance  from  the  geocenter  to  the  vehicle  CG,  measured  along  the  vector  whose 
orientation  is  determined  by  the  angles  8  and  X£  is 

R  =  R£  +  ALT  (3) 

where  R£  is  defined  by  equation  (1). 

It  may  be  worthwhile  to  digress  briefly  on  the  subject  of  latitude. 

There  are  three  latitudes  commonly  used  in  astrodynamics;  however,  in  this 

report  we  will  be  concerned  only  with  geocentric  and  geodetic.  Geocentric 

latitude,  X^,  is  defined  as  the  angle,  measured  at  the  earth's  center,  between 

the  equatorial  plane  and  a  radius,  R^,  to  the  point  of  interest  on  the  earth’s 

surface.  The  geodetic  latitude,  X  ,  is  the  angle  between  the  equatorial  plane 

S 

and  a  normal  to  a  reference  spheroidal  surface,  ZR,  as  shown  in  figure  3. 

The  reference  spheroid  is  an  oblate  ellipsoid  whose  surface  nearly 
approximates  the  sea-level  surface  of  the  earth.  Such  a  spheroid  is  character¬ 
ized  by  a  "flattening"  or  "oblateness,"  f,  where 

f  =  (a  -  b)/a  (A) 


and  a  is  the  semimajor,  and  b  the  semiminor  axis  of  the  ellipse  of  revolution 
forming  the  spheroid.  The  derivation  of  equation  (1)  can  be  easily  obtained 
from  the  equation  for  an  ellipse  if  a  equals  equatorial  radius,  R£q,  and  b 
equals  polar  radius,  RpQ.  Therefore, 


(5) 
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Solving  for  R^. 


(8) 


From  reference  10 


and  from  equation  (4) 

*Po  1-f 


(9) 


By  squaring  both  sides  and  substituting  equation  (9)  for  f. 


(10) 


Now,  by  adding  -1  to  both  sides. 


From  equation  (2) ,  it  is  shown  that 


P 


E 


(11) 


Equation  (1)  is  now  obtained  by  substituting  equation  (11)  into  equation  (8) 


R 


E 


EQ 


(l  +  Pr  sin2 
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Che  sea-level  radius  of  the  earth  can  now  be  determined  at  any 
station  by  ’'..uwing  the  geocentric  latitude,  X^,  of  the  station. 

The  astronomical  latitude,  X  ,  is  defined  as  the  angle  between 

3 

the  local  astronomical  zenith,  Z  ,  or  vertical  (as  determined  by  a  plumb  bob 
or,  actually,  by  the  direction  of  the  local  gravitational  field)  and  the 
equatorial  plane.  Since  astronomical  latitude  is  based  upon  the  local  gravita 
tional  fieli.'  and  is  also  centrifugal-force  affected,  it  differs  from  the 
geodetic  latitude  by  a  very  small  value  and  can  be  considered  negligible. 
Therefore,  this  report  will  consider  geodetic  and  astronomical  latitudes 
identical  and  will  be  concerned  only  with  the  relationship  between  geodetic 
ar.d  geocentric. 

To  convert  from  geodetic  latitude  to  geocentric  latitude,  the 
slope  (0)  of  the  tangent  (geodetic  horizon,  see  figure  4)  at  the  station  on 
the  earth's  surface  can  be  obtained  by  differentiating  equation  (5). 

2x  +  2z  dz  _  q 
d  2  R_  2  dx  " 

fcQ  rO 


Figure  4.  Geodetic  to  Geocentric 
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With  the  foregoing  relationships  established,  the  second 
reference  frame  can  now  be  defined. 

(2)  Topocentric  or  Horizontal  Reference  Frame 

The  second  of  these  coordinate  systems  is  consnonly  known  as  the 
topocentric,  but  has  other  names  depending  upon  the  users  application.  Some 
of  the  names  are  horizon,  launch,  reentry,  and  trajectory  system.  This  report 
refers  to  this  coordinate  system  as  the  horizontal,  H.  The  word  topocentric 
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is  derived  from  the  Greek  ix>cos,  meaning  a  place.  In  thi":  reference  frame, 
the  origin  is  taken  at  the  observer  as  shown  in  figure  5.  The  Y^-axis  is 
directed  toward  the  east,  the  X^-axis  toward  the  south,  and  the  Z^-axis  is 
directed  toward  the  astronomical  zenith.  This  system  is  rotated  along  with 
the  earth's  reference  frame,  and  is,  therefore,  not  inertial.  The  angle  from 
the  north  direction  measured  clockwise  (or  east)  in  the  X^,  plane  is  defined 
as  azimuth,  AZ. 

With  the  horizontal  reference  frame  now  defined,  it  is  now 
possible  to  complete  three  of  the  six-degrees-of-freedom. 

b.  Equation  Development 

Newton's  second  law  for  a  rotating  coordinate  system  is  repeated  here 
for  convenience. 


or 


m 


’T- 

d2R 

dt2 


/  ,  F  -  2mm  x  v  -  mm  x  (u  x  R) 


d2R 

dt2 


Zf 


(15) 


where 


the  vector  acceleration  of  the  vehicle  CG  with  respect  to 
the  inertial  reference  frame 

the  total  force  vector  on  the  reentry  vehicle 


the  coriolis  acceleration  of  the  axis  system  due  to  the 
rotation  of  the  earth 

vehicle  mass 

the  centrifugal  acceleration  of  the  axis  system  due  to 
the  earth's  rotation 


ti)£  =  0.7292  x  lG-i*  radian/sec,  earth  rotational  angular  velocity 


Zj.j  will  rotate  with 
the  earth,  several  vector  cross-products  will  appear  in  the  final  equation  for 
d2R/dt2.  Therefore,  the  rotation  rate  vector,  o>E,  and  general  vector,  R,  must 
be  expressed  in  the  same  reference  frame  Y^.,  Z^.j  for  expansion. 

12 


Since  the  horizontal  reference  system  ^X^,  Y^, 
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As  shown  in  figure  5,  pseudo-trajectory  reference  frame  (Xg,  Yg,  Zgj 
has  its  origin  at  the  initial  position  of  the  vehicle  center  of  gravity.  This 
frame  is  fixed  in  space  (with  respect  to  the  horizontal  (H)  frame),  with  Zg 
equal  to  the  reentry  altitude,  and  with  Xg  and  Yg  rotated  from  Xg,  Yg  about  Zg 
seme  initial  azimuth  angle,  AZ.  The  vehicle  flight  path  is  initially  fixed  in 
the  Y*  -  Z*  plane  directed  downward  some  reentry  angle,  y  . 

a  n  C 

Therefore,  the  components  of  Xg 


X*  =  cos  AZ  *  Xu  +  sin  AZ 
n  n 

•  Y 

H 

(16a) 

Yp  *=  sin  AZ  •  Xg  -  cos  AZ 

*  Y 

XH 

(16b) 

ZK=ZH 

(16c) 

As  previously  shown,  the  horizontal  reference  frame  is  not  inertial, 

— ► 

and  therefore,  rotates  at  the  same  rotation  rate  as  the  earth.  So  that  nay 
be  written 


»  ZH  °n  ^H*  ^*H’  ZH 


ur  =  u_  sin  X  Z.,  -  cos  X  Y„ 
EE  g  tl  E  g  n 


(17) 


Equation  (17)  may  be  rewritten  after  substituting  equation  (16b)  into  equation 
(17). 


sin  X  Zg  -  (i!r  cos  X  •  sin  AZ  Xy 


+  w  cos  X  •  cos  AZ  Y„ 
Eg  H 


(18) 


Therefore,  the  components  of  the  earth's  rotation  rate  in  the  horizontal 
reference  frame  may  be  written 


“^II" 

““e 

cos 

*  sin  AZ 

(19a) 

“g 

Weyh  = 

cos 

A 

•  cos  AZ 

(19b) 

g 

_ 

sin 

(19c) 

““h  " 

~E 

g 
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To  detenaine  tbe  general  vector,  R,  and  displacenent  accelerations 
the  vehicle  CG  Best  be  considered  displaced  fron  the  horizontal  reference 
frane  by  the  vector,  R,  as  shewn  in  figure  6. 


\ 


VEHICLE 


I 
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From  figure  6,  the  following  relationship  can  be  established: 


R  =  R  +  r 
eo 


r  can  be  written  in  terms  cf  its  components  on  the  horizontal  reference  frame 


?  -  X  +  Y  +  Z  ZH 


where  X,  Y,  Z  are  the  scalar  vehicle  CG  displacement  on  the  horizontal 
reference  frame.  R  can  now  be  written  as 


t  «  X  it*  +  Y  +  Z  2U  +  & 

H  H  H  eo 


Although  R  can  be  written  in  any  'reference  frame,  its  vector  components  must 
be  in  the  same  reference  frame.  Therefore,  the  equation 


$  =  R  Z 

eo  eo  c 


must  be  written  in  terms  cf  the  horizontal  reference  system. 

From  figure  4,  the  vector,  ?c,  in  the  geocentric  reference  frame  can 
be  written  in  the  form 


2  =  cos  n  +  Y„  cos  n 
c  H  H 


or  by  substituting  in  equations  (16b)  and  (16c),  the  following  is  obtained: 


Z  =  Zu  cos  n  +  XA  sin  AZ  sin  n  -  Y'  cos  AZ  sin  n 
C  n  rt  ri 


Equation  (23)  can  now  be  rewritten  as 


5  =  R  sin  AZ  sin  n  X'  -  R  cos  AZ  sin  n  Y'  +  R  cos  o  Z  (25) 

eo  eo  H  eo  H  eo  H 


where  components  of  F  in  horizontal  reference  frame  are 
r  eo 

R^  =  R  sin  AZ  sin  r, 

R,„,  =  -R  cos  AZ  sin  n 
YH  eo 

R,„  =  R  v.os  n 
ZH  eo 
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By  substituting  and  collecting  terms,  the  equation  for  K  is 

R  ■  (RxH  +  X)  +  (RyH  +  V)  %  +  (RzH  +  z)  2H  «' 

The  vectors  R  and  have  now  been  expressed  in  the  horizontal 
reference  frame.  The  development  of  equation  (15)  can  proceed  further  by 
taking  the  first  derivative  of  R,  with  respect  to  time,  in  the  following 
general  form  for  an  inertial  coordinate  system: 

V't  /  l_atj 


where 


is  the  derivative  relative  to  the  rotating  reference  system 


(d\ 

\dt  J 


is  the  vehicle  CG  velocity  vector  with  respect  to  the  inertial 
reference  frame  (R) 


Therefore, 


I  -  (RxH  +  R)  *H  +  (RyH  +  i)%*  (*zH  +  Z)  ZH 

+  *  (RxH  +  5)  ^  +  ( V  +  * )  ^  +  (RzH  +  Z)  <28) 

where  a  dot  above  the  variable  denotes  first  derivative  with  respect  to  time. 
The  last  term  may  be  expanded  in  the  following  matrix  form: 


X ' 
H 


w.  x  R  =  so. 
E  E 


E  E 

vH  zH 


j  (RxH  +  X)  (V  *  *)  (RzH  +  Z)  I 

■  “EvH  (RzH  +  2)  -  “E2„  ("yli  +  V)  SH+  -E2„(RxB  +  X) 
“  J  u 

-  -EkH  (RzH  +  z)]  %  +  [\n  (SVH  +  ')  -  -Evli  (R*,i  *  X) 


» 
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From  this  expanded  form,  the  following  equations  can  be  established: 


A  -  ^  (•* +  z)  -  x„  (V +  *) 


2H 


"  -  “Ez11  (RxH  +  *)  -  UExH  (*zH  +  Z) 


c  = 


\„(V  +  Y)"%h^  +  X) 


Hierefore,  the  velocity  vector  is 


=  (X  +  A)  X^  +  (Y  +  B)  Y^  +  &  +  C)  ZH 


To  obtain  the  acceleration,  the  second  derivative  of  R  with  respect  to 
taken 


d2R 


dt 


2  =  (X  +  A)  X^  +  (Y  +  B)  Y^  +  (Z  +  C)  ZH 


+  ajF  X  (X  +  A)  +  (Y  +  B)  +  (Z  +  C)  Z 


‘H 


H 


Sh] 


•  « 


where  A,  B,  C  are  the  derivatives  of  equations  (29a),  (29b),  and  (29c), 


A  =  io_  Z  -  tu  Y 
£  E 

yH  zH 


B  =  w_  X  -  w_  Z 
EzH  ExK 


C  =  oj_  Y  -  w_  X 
EXH  fcyH 


and  the  last  term  is  expanded  to  give 


(29a) 

(29b) 

(29c) 

(30) 

time  is 

(31) 

(32a) 

(32b) 

(32c) 
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-  Y  dR 
“EX  dF 


XH  YH  ZH 


“e  we  we 

xH  yH  zH 


(X  +  A)  (Y  +  B)  (Z  +  C) 


=  U_  (Z  +  C)  -  10-  (Y  +  B)  Xl 
vH  EzH  _  11 


Kh  *  + 


A)  -  u  (Z  +  C)  Y£ 
xH  J 


f“E  „  «  +  B)  '  “E  „  «  +  A)1  % 
[  xH  yH  J 


By  substituting  the  above  into  equation  (31),  the  following  is 
obtained  for  the  acceleration 


[X  +  w£  Z  -  w£  Y  +  w£  (Z  +  C)  -  uE  (Y  +  B)1  X^ 
yH  zH  yH  zH  j 


+  Y  +  w£  X  -  u£  Z  +  w£  (X  +  A)  -  w£  (Z  +  C)  I  Y* 
L  zH  xH  zH  xH  J 


:z  +  C)  j  Y^J 


[Z  +  u£  Y  -  U  X  +  UJ£  (Y  +  B)  -  «E  (X  +  A)j  Z^  (33) 
xH  yH  xH  yH  J 


By  collecting  terms  and  reducing  equation  (33),  the  final  development 
of  vehicle  acceleration  with  respect  to  the  horizontal  reference  frame  is 


.(;  +  m  +  %hc-%hb)x. 
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jY  +  ?.B  +  uc  A  -  u  C)  Y’ 
\  “zH  xH  / 

(  Z  +  2C  +  (o_  3  -  A  \  Z„ 

V  ExH  EvH  /  H 


or  in  component  form  in  the  horizontal  reference  frame. 


a  =  X  +  2A  +  u)_  C  -  to  B 
EyH  EzH 


(35a) 


a  =  Y  +  2B  +  io_  A  -  to_  C 
y  zH  xH 


(35b) 


Z  +  2C  +  to_  B  -  to_  A 
xH  EyH 


(35c) 


This  completes  the  derivation  of  the  acceleration  equations,  which 
includes  coriolis  and  centripetal  acceleration  due  to  rotation  of  the  coordinate 
system,  and  accelerations  of  the  vehicle  CG  due  to  forces  and  moments. 

To  complete  the  equations  of  motion,  it  will  be  necessary  to  introduce 
the  forces  and  moments  acting  on  the  vehicle.  These  forces  and  moments  result 
mainly,  for  this  program,  from  gravitational  attraction,  and  from  reaction  of 
the  air  on  the  vehicle  as  a  result  of  its  motion.  This  is  not  to  say  other 
types  of  forces  cannot  be  involved. 

c.  Earth  Position  Coordinates 

Since  the  components  of  the  gravitational  force  are  dependent  upon  the 
position  of  the  vehicle  in  the  earth  coordinate  system,  the  instantaneous 
longitude  and  geocentric  latitude  can  be  obtained  by  the  following  method. 

From  equation  (6)  and  reference  9,  the  geocentric  latitude  is 


X  ~  arcs in 
c 


and  the  longitude  is 
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6  =  arctan  — 
XE 


where  X^,  Y^,  and  Z,,  are  coordinates  in  the  earth  reference  system. 

This  program  expresses  the  vehicle  position  in  the  horizontal  reference 
frame;  therefore,  the  position  must  be  transformed  to  the  earth  reference  frame. 
This  can  be  done  through  examination  of  figures  2,  5,  and  6,  and  the  following 
transformation  matrixes. 


to 


y 

“E* 


The  order  of  transformation  will  be  from  X’, 

H 

V  nnA  7 

*E* - E* 


ZH  CO  XH>  YH’ 


and  Z 


H’ 


Therefore, 


Xu  =  X*  sin  AZ  -  Y'  cos  AZ 

n  ri  n 

Yr  =  X^  cos  AZ  +  Y^  sin  AZ 
ZH  "H 

or  in  matrix  form 


1 

X 

re 

sin  AZ 

(-cos  AZ) 

o 

XE 

yh 

= 

cos  AZ 

sin  AZ 

0 

• 

yh 

ZH 

0 

0 

1 

ZH 

(36a) 


(36b) 


and 


X  =  X„  sin  X  +  Z„  cos  X 
E  li  g  H  g 


Y  =  Y 

E  il 


Z  *  -X,,  cos  X  +  Z.,  sin  Xg 
E  II  g  ri 


(37a) 


from  which  the  following  matrix  is  formed: 


21 
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X  sin  X  0  cos  A 

e  g  g 

Y  0  10 


Z  -cos  X  0  sin  X 

,  ej  L  g  8. 


Fxh] 


(37b) 


By  substituting  equation  (36b)  into  equation  (37b) 


sin  X  0  cos  X 
g  8 


0  1 


sin  AZ  (-cos  AZ)  0 


cos  AZ  sin  AZ  0 


Z  -cos  X  0 

L  e_J  L  8 


sin  X  I  I  0 


and  expanding  equation  (38), 


(sin  Ag  sin  A.z)  (-sin  Ag  cos  Az)  cos  X? 


cos  AZ 


sin  AZ 


(-cos  Xg  sin  AZ^  (cos  Ag  cos  AZ^  cos  Xg  Z^ 


From  equation  (39) 


=  X *  (sin  Ag  sin  Az)  -  Y^  (sin  Ag  cos  Az)  +  Z’ 


•  COS  A 

H  g 


(40a) 


where 


Y_  =  X*  cos  AZ  +  Y’  sin  AZ 
t  H  n 


XH  =  RxH  +  X 


YH  RyH  +  Y 


(40b) 


=  -X'H  (cos  Ag  sin  AZ^  +  Y^  (cos  Ag  cos  Az)  +  Z*  sin  Xg  (40c) 


Z’  =  R  „  +  Z 
H  zH 
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RxHf  RyH*  and  RzH  are  defined  by  equation  (26).  Therefore,  the  instantaneous 
geocentric  latitude  can  be  obtained  from 


X 

c 


(41) 


where 


R  = 


1/2 


since  the  earth  is  assuned  to  be  a  body  of  revolution  about  the  polar  axis,  the 
instantaneous  longitude  may  be  defined  as 


6  « 


(42) 


and  will  be  considered  positive  when  measured  east  from  the  prime  meridian  and 
negative  when  measured  west. 

d.  Gravitational  Equations 

The  potential  function  and  equations  used  to  obtain  the  gravitational 
forces  to  be  included  in  the  equations  of  motion  are  presented  in  this  sub¬ 
section.  The  gravitational  potential  of  the  earth  is  an  expansion  of  spherical 
harmonics  as  a  function  of  latitude  and  longitude.  However,  in  this  program, 
the  earth  is  assumed  to  be  a  body  of  revolution  so  that  the  longitudinal  effects 
may  be  ignored.  Although  some  satellite  data  indicate  that  the  earth’s 
equator  is  elliptical  where  the  major  axis  is  somewhat  larger  than  the  minor 
axis,  the  difference  though  is  very  small  compared  to  the  equatorial  radius. 

The  gravity  components  are  determined  by  taking  partial  derivatives  of 
the  gravitational  potential  equation  (Ref.  9'. 


(43) 
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where 

G  M  =  1.4076427  x  10lfc  ft3/sec2  gravitational  constant  of  the 

O  O  . 

earth 

R  =  radial  distance  to  body  CG  from  earth’s  geocenter 
R_  =  earth's  equatorial  radius 
Xc  =  geocentric  latitude 

J  =  1623-41  x  10_fe,  coefficient  of  the  second  harmonic 
H  =  6-04  x  10-k ,  coefficient  of  the  third  harmonic 
K  =  6.37  x  10“°,  coefficient  of  the  fourth  harmonic 


The  gravitational  potential  field  is  usually  defined  by  spherical 
harmonics.  Each  higher  order  term  in  equation  (43)  is  due  to  a  deviation  of 
the  potential  from  that  of  a  true  sphere-  Therefore,  the  potential  is  expressed 


as  a  series  of  Legendre  polynomials  in  the  variable,  X  .  Many  analyses  (Refs. 

9,  11,  12)  consider  only  the  second-,  third-,  and  fourth-order  terms.  The  first 


harmonic  is  missing  because  of  the  choice  of  an  equatorial  coordinate 

system  in  which  X  is  measured  relative  to  the  center  of  the  earth.  The  fifth 


harmonic  is  not  sufficiently  well  known  to  justify  its  inclusion  in  this 


equation. 


The  gravitational  acceleration  of  a  vehicle  along  any  line  is  the 
partial  derivative  of  the  potential  function,  v,  along  that  line.  See  figure  7. 


Fieure  7.  Gravitational  Direction 
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g 


is  related  to 


s  -  -«(*.  0 

where  7$  is  the  gradient  vector  of  v^R,  X  from  figure  7,  the  force  along 
the  radius,  R,  is 


G  =  mg  =  -m 
r  ®r 


and  in  R,  X  direction 
c 


c 


where  n  =*  mass  of  vehicle  in  slugs.  From  equation  (43),  G  and  G*  are 

K  C 

obtained. 


simplifying  and  dividing  through  by  c. 


23 


AFWL-TR-68-61 


+  *1  fcfi)'  (3  sin  X  -  5 

3  \  R  /  '  c 


1-3  sin2  A 


2>c) 


5  sin^  X 


(-p)  (3  -  30  sin2  Ac  4 

“■f-'Wi- 


+  35  sin4  A 


G  M 

GA - 2_o  _2 

C  R2 

L 


sin  A  cos  A 
c  c, 


(44a) 


(cos  A£  -  5  sin2  A^  cos  A£^ 

(-3  sin  A^  cos  A£  +  7  sin3  A£  cos 


(44b) 


Since  G  is  along  the  radius  ^  and  GA  is  perpendicular  to  R,  they  are  in  the 
inertial  reference  frame.  For  them  to  be  included  in  the  equations  of  notion, 
their  components  must  be  found  with  respect  to  the  horizontal  reference  frame. 
This  can  be  accomplished  by  the  use  of  direction  cosines  and  transformation 
matrices. 

After  examination  of  figure  8,  the  following  series  of  equations  can 

be  written  to  transpose  GD  and  GA  to  the  horizontal  reference  frame. 

K  C 

It  is  noted  that  G  and  GA  are  taken  as  positive  in  the  direction 

K 

shown  in  figures  7  and  8.  Therefore,  by  definition 


GR=+5R 


1  34 
G>  "  +  R  3R 


G5  =  0 
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Equations  (46)  and  (47)  can  be  expressed  in  matrix  form 


GX, 


GY, 


gzk 


(-sin  X  cos  X  cos  6  +  cos  X  sin  X  ) 
\  g  c  g  c/ 


^-cos  X£  sin  gj 


(-cos  X  cos  X  cos  8  -  sin  X  sin  X  ) 
\  g  c  g  c/ 


\ 

/  \ 

r  i 

(sin  X 

sin 

A  cos 

8  + 

cos 

X 

cos  X  | 

[sin  X  sin  gj 

GR 

V  g 

c 

g 

cy 

\  g  / 

1 

^sin  Xc 

sin 

B) 

(-cos  8) 

• 

GX 

/ 

jcos  X 

sin 

X  cos 

8  - 

sin 

X 

cos  X  i 

(cos  X  sin  8l 

G8 

V  g 

c 

g 

c) 

\  g  /  J 

(49 


2  written  into 

the 

following  equations  for 

Gil  = 

-sin 

X 

CCS 

X 

cos 

8  + 

cos 

X 

sin 

X 

g 

c 

g 

G12  - 

sin 

X 

sin 

X 

cos 

8  + 

cos 

X 

cos 

X 

g 

c 

g 

G21  = 

-cos 

X 

c 

sin 

6 

G22  - 

sin 

X 

sin 

8 

c 

G31  - 

-cos 

X 

cos 

X 

cos 

8  - 

sin 

X 

sin 

X 

g 

c 

g 

G32  = 

+COS 

X 

sin 

X 

cos 

6  - 

sin 

X 

cos 

X 

g 

c 

g 

t 

(so: 


The  final  cransformation  is  obtained  by  substituting  equation  (49) 
into  equation  (48) . 


GX^  =  GR  (Gil  •  sin  AZ  +  G21  •  cos  AZ) 


+  GX  (G12  •  sin  AZ  +  G22  •  cos  AZ) 


GY’  »=  GR  (-G11  •  cos  AZ  +  G21  •  sin  AZ) 

n 


+  G'  (-G12  -  cos  AZ  +  G22  •  sin  AZ) 


GZ'  =■  GZU  =  GR  (G3i)  +  GX  (G32) 
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The  components  of  the  gravitational  forces  are  now  expressed  in  the 
horizontal  reference  frame  and  may  be  included  with  equations  (35)  to  complete 
the  equations  of  motion. 

Before  the  aerodynamic  forces  and  moments  can  be  resolved,  a  method 
must  be  established  for  the  transformation  of  components  from  the  body 
reference  frame  to  the  horizontal  reference  frame  and  the  inverse. 

e.  Euler  Angles 

Up  to  this  yoint,  only  the  equations  for  the  three-translational 
degrees-of-freedom  have  been  derived.  This  establishes  the  vehicle  CG  with 
respect  to  a  horizontal  reference  system.  To  analyze  the  body  forces  and 
moments  as  a  result  of  its  motion  about  the  CG,  the  angular  orientation  of 
the  vehicle  must  be  known  with  respect  to  the  horizontal  reference  system. 

The  Euler  angles  describing  the  angular  motion  are  <j>,  6,  and  y,  and  will 
be  noted  as  pitch,  yaw,  and  roll  angles,  respectively.  These  angles  are  to  be 
differentiated  from  the  angles  of  attack  and  sideslip.  The  angles  of  attack 
and  sideslip  have  a  kinematic  definition  based  on  the  components  of  the  free- 
stream  velocity  relative  to  the  body  axes  of  the  vehicle.  The  angular  displace¬ 
ments,  6,  and  i|>,  on  the  other  hand,  are  used  to  measure  the  vehicle  attitude 
with  respect  to  a  fixed  set  of  axes  (X,  Y,  and  Z) ,  and  in  no  way  require  motion 
of  the  vehicle  relative  to  the  surrounding  air  for  their  definition. 

Consider  a  vehicle  moving  with  respect  to  the  inertial  axes  X,  Y,  and 
Z  which  are  also  fixed  ir.  space.  The  following  will  describe  one  of  several 
ways  of  specifying  the  angular  orientation  of  a  vehicle  at  any  instant  of  time. 
This  winl  be  done  by  successively  pitching,  yawing,  and  rolling  the  X,  Y,  and 
Z  axes  until  they  coincide  with  the  axes  x,  y,  and  z  of  the  vehicle  as  shown 
in  figure  'J. 

First,  pitch  the  vehicle  by  an  angular  displacement  <$  ,->round  OX  so  that 
Z  goes  to  Zj  and  Y  goes  to  Yj  (figure  9a).  The  relationship  between  the  two 
coordinate  systems  is  then  given  by  the  following  transfer  matrix: 
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Next,  yaw  the  vehicle  by  an  angle  0  about  the  OYj  axes  so  that  Xj  goes 
to  X2  and  Zj  to  Z2  (figure  9b).  The  new  relationship  is  given  by  the  transfer 
matrix 


cos  6  0  -sin  0 

0  10 


Z2  I  sin  0  0  cos  0  j  Zj  j 

»-  —  l—  —  _ ■ 


M 

i 


Finally,  roll  the  vehicle  by  angle  around  the  OZ?  axis  so  that  X2 
moves  to  X3,  and  Y2  to  Y3  (figure  9c).  The  transfer  matrix  for  this  rotation 
is 


H 


cos  tj;  sin  0 

-sin  £  cos  0 

0  0  1 


where 


x  =  X. 


y  -  y. 


z  =  Z. 


The  operations  of  pitch,  yaw,  and  roll  for  th.'.s  program  must  always 
be  performed  in  this  order  since  angular  displacements  do  not  follow  the 
ordinary  laws  of  vector  addition,  but,  in  fact,  follow  a  nonccmmutative  law. 
Under  this  system  of  rotations,  the  direction  cosines  of  the  finial  vehicle 
axes  x,  y,  and  z  to  the  fixed  axis  X,  Y,  and  Z  are  obtained  by  substituting 
equation  (52)  into  equations  (53). 


1 — 

<N 

X 

1 _ 

cos  6 

(-sin  0 

sin  ^) 

(-sin 

cos  ; )  | 

!  X  j 

i  ! 

?2 

= 

0 

'■'OS 

- 

-sin 

a  i 

| 

*  i Y ! 

Z= 

L  — 

sin  0 

(sin  t 

COS  0) 

(cos  o 

COS  ••) 

1  z 

1  1 

(55) 
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Then  substitute  equation  (54)  into  equation  (55) 


(cos  ip  cos  0)  (-sin  0  sin  $  cos  ip  +  sin  ip  cos  6,» 
(-sin  <p  cos  9)  (sin  0  sin  $  sin  £  +  cos  cos  <p) 

sin  0  (cos  0  sin  <p) 


(-sin  9  cos  <)>  cos  y  -  sin  y  sin  i>) 
(sin  ij;  sin  9  cos  <£  -  sin  $  cos  ij>) 
(cos  0  cos  <p) 


1 


(56) 


Therefore,  equation  (56)  will  transform  components  from  the  fixed  axes  or 
horizontal  reference  system  to  the  vehicle  or  body  reference  system.  To  trans¬ 
form  components  from  the  vehicle  axes  to  the  horizontal  references  axes,  the 
matrix  of  equation  (56)  must  be  inverted  by  interchanging  rows  and  columns. 

The  elements  of  equation  (56)  may  be  written  for  simplicity  as  follows,  where 
subscripts  indicate  rows  and  columns. 

Aj ,  «■  cos  0  cos  p 

A12  «  -sin  0  sin  $  cos  <p  +  sin  y  cos  6 

Aj.  *=  -sin  0  cos  cos  $  -  sin  ip  sin  $ 

A2  j  *  -sin  ip  cos  0 

A22  *  sin  0  sin  0  sin  ip  +  cos  ^  cos  $ 

J,.  =  sin  ifi  sin  9  cos  $  -  sin  $  cos  <P 

sin  & 
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cos  0  sin  6 
cos  0  cos  <P 


(57) 


f.  Transformation  of  Angular  Velocities 

It  is  now  necessary  to  express  the  angular  velocities  w  »  and  w. 
with  respect  to  the  vehicle  axes  in  terms  of  the  Euler  angles. 
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2 


*2 


Figure  10.  Angular  Velocities  of  Euler  Angles 


The  components 
to  the  vehicle  axes  x. 


of  the  angular  velocities  along  the  X,,  Y?, 
y,  and  z  shewn  in  figure  10  are 


and  Z,  axes 


u 

x 


y 


dS 

dt 


cos 


sm 


.  dc 

+  dFcos 


-dt 

dt 


sin 


+ 


do 

dt 


(58) 


The  equations  in  section  Il-2e  and  equation  (56)  have  been  derived 
with  respect  to  the  inertial  axes  X,  Y,  and  Z.  If  the  axes  system  is  allowed 
to  rotate,  the  angular  velocity  of  the  earth  must  be  included.  Therefore, 
the  inertial  angular  rates  of  the  vehicle  are 
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“Y  +  "E.. 


n2  =  WZ  +  UE„ 


where  w  ,  w  ,  and  u_  are  defined  by  equation  (19).  These  are  the 
xH  yH  zH 

components  of  the  earth’s  rotation  rate  in  the  horizontal  reference  frame. 

Since  this  program  was  written  for  a  symmetrical  cone  with  uniform 
weight  distribution  and  further  assumes  symmetrical  ablation,  tnen  the  products 
of  inertia  will  be  zero  regardless  of  roll  angle,  si.  (This  will  be  shown  ir« 
section  II-2g.)  Therefore,  the  roll  and  transverse  moment  of  inertia  can  be 
obtained  independent  of  This  will  simplify  the  equations  and  allow  the 
angular  momentum,  hQ,  to  be  written  in  the  second  rotation. 

Therefore,  equation  (58)  can  be  simplified  and  written  in  the  X9,  Y2, 
and  Z2  axes  as 


dO  n 

03  =  -  —  COS  9 

X2  dt 


V2  dt 


WZ2  =  '  H  Sin  6  *  If 


which  arc  tae  total  rotation  rates  of  the  vehicle  in  the  second  rotated  axes. 

Now  03-  ,  u>_  ,  and  a-  must  also  be  written  in  the  same  axes  system 

ScH  xH  zH 

as  the  vehicle  angular  rates.  This  can  be  done  by  using  the  transformation 
equation  (55),  which  will  result  in 


“E  E 


cos  8  -  oiv  sin  6  sin  $  -  u  sin  0  cos  c 


=  0  +  03-  cos  o  -  w_  sin 
EE  E_ 


03-  =  os-  sin  6  +  sin  c  cos  0  +  cos  c  cos 

t.  t  E  E 
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Therefore,  equation  (59)  can  be  written  as 


f;  =  u  + 

x  x2  Ex, 


£1  =  u  + 

y  y2  £y2 


—  Lj  +  U- 
z  Z2  Ez2 


which  are  the  inertial  angular  velocities  of  the  vehicle  in  the  X,,  Y2,  and  Z, 
axes. 

Since  equation  (62)  is  written  in  the  second  rotated  axes  system,  it  is 
also  necessary  to  have  the  angular  velocity  of  that  coordinate  system.  There¬ 
fore, 

H  cos  6  +  io 

x.  dt  E.. 

/  A, 

n  =  |§-  +  u 
y2  dt  eY2 

£5  =  -  4f  sin  6  +  w  (63) 

z2  dt  Ez 

2 

which  completes  the  transformation  of  angular  velocities, 
g.  Moment  of  Momentum 

A  particle  of  mass,  m,  moving  with  velocity,  R,  figure  11,  has  a  linear 
momentum  (Ref.  13) 

-*• 

P  =  a  R 

The  moment  of  this  linear  momentum  about  an  arbitrary  point,  0,  is 


defined  as 


where 


h  =  r  x  m  R 
o 


R  =  absolute  velocity  of  m 

r  =  radius  drawn  from  0  as  shown  in  figure  11 
hQ  =  the  angular  momentum  of  the  particle  about  the  point  0 
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Figure  11.  Moment  of  Momentum  about  0 

The  angular  momentum  of  a  rigid  body,  figure  12,  about  any  axis  perpen¬ 
dicular  to  the  plane  of  motion  and  passing  through  the  point  0  fixed  in  a 
moving  body  is  the  sum  of  the  moments  of  the  linear  momenta  of  ail  its  particles 
abcut  the  axis.  Therefore,  equation  (64)  can  be  rewritten  as 

h  =  >  r ,  x  m  v, 

o  i  i 

where  =  velocity  of  any  point  i  on  the  body.  The  velocity  of  a  representa¬ 
tive  particle  of  mass  ck  may  be  expressed  in  terms  of  the  velocity  of  0  plus 
the  velocity  of  m^  with  respect  to  0. 

v.  =  v  +  w  x  r . 
i  o  l 


Therefore, 


h 

c 


x  r 


*) 


or  is  expanded  from 
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Figure  12.  Coaponents  of  Moaentua 


If,  by  definition,  0  coincides  with  the  center  of  mass,  the  tera  a.  r.  is 
equal  to  zero,  the  angular  aoaentin  can  be  expressed  by  the  following  integral 

h  *>  1  r  x  (w  x  r)  da  (65) 

°  J 

The  relation  between  the  angular  aoaentua  of  a  body  and  the  applied 

-•  •• 

aoaents,  M  ,  of  the  forces  /  .  F  »  aR  acting  on  a  is  obtained  froa  rotational 
equation  of  aotion. 


M  «  T*  r  x  aR 


(66) 


f 
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Taking  Che  second  derivative  of  the  expression 


R  =  R  +  r 
o 


and  substituting  into  equation  (66) 


Mo  =  ^  ^  x  m  (rq  +  r) 
or 

**o  =  £  (r  x  m  r)  -  Rq  x  £  nr  (67) 

Again  mr  *  0  when  0  coincides  with  center  of  mass.  Equation  (67)  reduces 
to 


M 

o 


£ 


d 

dt 


-*■ 

(r  x  m 


r) 


(68) 


but  r  is  equal  to  the  following: 


where 


total  derivative  relative  to  inertial  axis 


derivative  relative  to  rotating  axis  and,  since  body  is 
rigid,  is  equal  to  zero 


There fore. 


r  “  x  r 


and  substituting  into  equation  (67)  the  relationship  between  M  and  h  is 

c  o 

established. 


V? 

o 


dt  V 


r  x  a  U  x 


(69) 


38 


AFWL-TR-68-61 


or 


M 

o 


d  h 

o 

dt 


(70) 


With  this  relationship  established,  the  derivation  of  hQ  can  be  continued.  Th> 
integral  in  equation  (65)  nust  be  evaluated  and  this  is  done  by  first  expanding 
the  cross  product  <u  x  r 


it  x  r 


j 


j 

k 


(71) 


where  and  H 2  is  defined  by  equation  (62).  Then  expansion  of  the  cross 


product  X  x  (ax  r)  is 


r  x  (a  x  r)  = 


I  (V  -  V)  (w)  (w) 


i  (y2  +  z2)  -  f.v  (xy)  -  (xz)j  d= 
3  [~nx  <*‘>  +  E„  ♦  z2)  ~  -z  (>*=)] 
k  (xz)  +  (yz)  +  (x-  +  y~)J 


cc 


(72) 


3v  definition,  the  nccents  of  inertia  of  a  body  about  the  x,  y,  z  axes  are 

i  -  (  (>*'  +  z:)  <=• 


TSTr 
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and  the  products  of  inertia  are 


V  ) 

xy 

dm 

1  =  f 

xz 

dm 

xz  j 

1  =  f 

yz 

dm 

yz  ] 

If  the  x,  y,  z  axes  are  considered  the  principal  axes,  the  products  of  inertia 
and  the  time  derivatives  of  the  moments  of  inertia  are  zero. 

Therefore,  by  integrating  equation  (72)  over  the  entire  body  and 
applying  the  above  definition,  the  following  is  obtained: 


h  i  +  Ifl  j  +  I  Q  k 

o  xx  y  y  J  zz 


(73) 


in  which  case  the  moment  of  momentum  components  along  the  x,  y,  z  axes  are 


h  =  I  ft 

X  XX 


h  *  I  ft 

y  y  y 


h  =  I  ft 
z  z  z 


(74) 


It  was  shown  previously  that  the  moment  about  a  body's  center  of  gravity  was 
equal  to  the  time  derivative  of  the  moment  of  momentum. 


d  h 

M  =  — - — - 
o  dt 


which  in  general  form  can  be  expressed  as 

/ dho\  _  dhol 

\dt  f  ~  [dt  J 

Substitute  this  into  equation  (70),  Euler's  moment  equation  is  obtained. 


+  u  x  h 


8  =|S*s,«h 

o  dt  2  o 


(7r>) 


where  ft  is  the  angular  velocity  of  the  second  rotated  axes  system. 
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Taking  the  cross  product  x  h  in  a  manner  similar  to  equation  (71), 


equation  (75)  becomes 


M  =  fe  i  +  d^J  +  |l£^l+/Q  h  -  ft  h  ) 
o  |_  dt  dt  J  dt  J  ^  y2  z  22  y) 

+  (°«h*  -  Wz)  5  +  (Wy  '  V*)  E 


and  collecting  terms  yields 

m  -fei+n  h  -c  h)i  +  fe*  +  c  h  -si  h)j 

o  \dt  y?  2  y2  y J  Vdt  22  x  X2  2/ 

+  +  fi  h  -  fi  h\k  (77) 

\dt  x2  y  y2  xy 

Now,  taking  the  derivatives  of  equation  (74)  and  remembering  d(I)/dt  =  0  for 
principal  aces, 


fa.ii 

dt  xx 


i  h 

dt  y  y 


?-ID 

dt  y  y 


where  are  the  derivatives  of  equation  (62) 


S?  =  cos  6  +  6  e  sin  9  +  w.. 
x  r,x 


H  =  $  -f  w 

y  -y2 


Q  =  sin  6  -  c  0  cos  B  +  c 
2 


and  to,.  ,  to—  ,  to_  are  the  derivatives  of  equation  (61).  Before  taking  the 
Ex-j  Ev~  Ez? 

derivatives  of  equation  (61),  let  the  following  be  elements  of  equation  (55) 
where  subscripts  indicate  row  and  columns: 


i 


(80) 


(81) 


(82) 
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Substituting  equations  (74'  (7£>),  and  (79)  into  equation  (77)  and  writing  Mc 

in  component  form  with  1^  =  1^, 

M  =  I  (— t  cos  6  +  J  6  sin  6  +  u.  j  +  H  „  (i  n  ^  -  Q  .  (in') 

X  x  l  v  Ex2  I  y2  \  z  z/  z2  V  x  y / 

M  =1  (  6  +  tl>_  )+n  (i  !!  )  -  !1  ( I  fi  \ 

y  x  l  Ey  /  Z2  \  x  x /  X2  \  z  z / 


M  =  I  (-£  sin  6  -  $  6  cos  6  +  !|>  +  w_  J  +  Q  (i  fi  ) 
Z  Z  \  Ez2^  X2  ^  X  yj 


-s,  (is  ) 

y2  \  x  x/ 


h.  Summary  of  Equations  of  Motion 

For  convenience,  the  equations  of  motion  are  restated  here  with  all 
terms  added  and  in  a  form  similar  to  that  used  in  the  program. 

=  £i-2A-a>  C  +  w  B-Gx 
dt2  m  EyH  EzH 


dk  =  IS.  _  28  -  A  +  C  -  Gy 
dt>  m  EzH  ExH 


d^z  Fz  _• 

-  - - 2C  -  ,  w  .  _ 

m  E„B+_  A-Gz 

dt-  xH  E  .. 

yH 


,2,  ..  n  i  n 

d^S  My  .  0  .  x2  z  z 

-  =  -=*-  -  9.  n  +  - - - u 

j .  2  1  Z2  x  I  Ey„ 

dt^  x  x  ^2 


(?* 


^4^  -  £  ,0  +  £  If  sin  e-  +  ; 
I  z~  y  dt  dt 

x 


d^o  Mz  "x2^x"v  ”y2^x”x  d‘ i 

-  =  - - — ■*-  -j-  — <■ — ; — —  4-  -  sm 

dt-  z  z  z  dt‘ 


.  di  dO 

T  IT  dT  cos  •-  “  -e 


rcyn  I1 . Mflilli'W'UWP1  J  .  !  !  W™ 5 .  ■ 11 


lAjjyjtm 


l 
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where  all  term  are  known  from  the  initialization  or  froir.  the  previous  time 
step. 

Although  these  equations  have  been  written  for  a  conical  shaped  body 
with  weight  change  due  only  to  ablation,  it  only  a  natter  of  modifying  the 
equations  to  include  thrust  and  additional  forces  and  moments  associated  with 
a  finned  missile. 

3.  INTEGRATION  TECHNIQUE 

To  integrate  a  set  of  n  simultaneous  second  order  differential  equations 
of  the  form 


*>’  v  •  •  •  <)■ 1  *  (h 

The  Runge-Kutta  method  of  numerical  integration  is  recommended  for  this 
type  of  application  and  especially  where  a  high  degree  of  accuracy  is  desired. 
Thus,  the  equations  of  motion  in  the  preceding  subsection  are  shown  as  six 
second-order  differential  equations,  which  are  to  be  solved  simultaneously. 

This  method  consists  of  four  sets  of  equations  which  involve  different  substi¬ 
tutions  into  the  differential  equations.  The  equations  are  solved,  and  the 
increments  of  the  functions  are  calculated  as  weighted  averages  of  the  solutions 
to  the  four  equations.  Therefore,  the  integration  procedure  for  the  expression 
of  the  form 


=  f  (t,  x,  x') 
dt2 

is  described  by  the  following  generalized  equations 


X  =  X  +  x' 
n+1  n  n 


XUl  =  Xn  +  r(k0  +  2ki  +  2k3  +  k2) 


where 


k  =  at  f  |t  ,  x  ,  x') 

0  V  r>  n  n/ 


k,  =  At  f 


ft  +  —  x  —  x’  x’  + 

V  n  2  »  Xn  ’  2  V  n  2  / 
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*  At  f  ft  +  At,  x  +  x*At  +  4^-  k, ,  x* 
\  n  n  n  2  1  *  n 


In  integrating  a  differential  equation  numerically,  the  equation  is  replaced 
by  a  difference  equation  and  solved  accordingly.  The  Runge-Kutta  method  has 
excellent  properties  for  stability  in  the  integration  and  if  the  integration 
step  is  taken  small  enough,  the  difference  equation  is  usually  close  to  the 
differential  equation  solution.  However,  some  variables  in  the  equations  of 
motion  may  have  peculiar  oscillations  or  increase  very  rapidly  which  could 
cause  an  unstable  condition  or  an  error  in  the  solution.  An  unstable  condition 
or  solution  can  often  be  made  stable  by  using  a  smaller  integration  step.  The 
reason  for  this  is  that  different  step  sizes  change  the  parameters  relating  to 
the  stability  of  the  difference  solution. 

The  integration  step  size.  At,  for  this  program  is  fixed,  although  the 
user  could  modify  the  integration  method  to  use  a  variable  step  size  under 
error  control  (Ref.  14).  The  user  has  the  option  to  break  into  the  program  at 
exact  specified  values  of  the  independent  variable,  t,  for  the  printout  step. 

4.  AERODYNAMIC  HEATING  AND  ABLATION 

The  equations  of  heat  transfer  are  general  in  that  they  apply  to  most 
axially  symmetric  vehicles  of  any  composite  skin  structure.  The  vehicle  skin 
is  assumed  to  be  made  of  discrete  layers  of  material  whose  properties  may  vary 
from  layer  to  layer  and  there  is  no  limit  on  the  number  of  layers.  This  program 
is  specific  in  that  the  equations  have  been  applied  to  a  blunt-nosed,  conical 
vehicle. 

Hypersonic  vehicles  may  be  categorized  into  two  general  types.  The  first 
consists  of  those  vehicles  which  must  be  propelled  through  the  atmosphere  to 
sustain  flight.  The  second  type  comprises  those  that  have  a  vast  store  of 
kinetic  and  potential  energy  which  will  be  dissipated  to  the  surrounding  air 
during  reentry.  The  latter  type  is  of  primary  interest  and  for  which  this 
program  was  written,  although  it  will  handle  either  type. 

In  either  case,  some  of  this  energy  is  expended  against  drag  forces  which 
assume  the  form  of  skin  friction  and  pressure  drag  which  produce  heat.  The 
heat  equivalent  of  the  kinetic  energy  contained  within  a  body  is  approximately 


t 
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V2  50,000  (Btu/lb) 
en 

where  V^is  the  reentry  velocity  in  ft/sec.  While  only  a  fraction  of  this 
energy  is  actually  transferred  to  the  body  as  heat,  this  fraction  is  still  a 
significant  quantity  and  accurately  illustrates  the  magnitude  of  the  problem. 

a.  Heat  Flux  Equation 

Experiments  in  high  speed  flow  have  verified  that  the  magnitude  and 
direction  of  heat  flow  at  the  surface  does  not  depend  on  the  difference  between 
the  wall  temperature  and  the  free-stream  temperature  as  in  low-speed  flow,  but 
rather  on  the  difference  between  the  wall  temperature,  TW,  and  the  adiabatic 
wall  temperature,  TAW.  It  Is  apparent  that  the  determination  of  the  adiabatic 
wall  temperature  will  be  of  prime  importance  in  the  calculation  of  heat  transfer, 
since  the  transference  of  heat  to  or  from  the  wall  will  depend  upon  whether  the 
skin  temperature  is  above  or  below  TAW.  The  adiabatic  wall  temperature  can 
conveniently  be  expressed  in  terms  of  a  dynamic-temperature  rise 


where 


(84) 


TEMPE  *  boundary  la  er  (B/L)  edge  temperature,  °R 

Vg  =  B/L  edge  velocity,  ft/a ec 

gc  =  gravitational  constant,  ft/sec2 

J  *  Joule’s  constant  =  778 

C  =  specific  heat  of  air,  Btu/lb-°R 
pa 

RF  =  recovery  factor  which  is  a  measure  of  the  fraction  of  the 
free-stream  dynamic  temperature  rise  recovered  at  ihe  wall 

Reference  15  shows  that  for  practical  purposes 


RF  =  yf~W- 

for  laminar  flow  and 


(85a) 


RF  nr  ^Tr 

for  turbulent  flow  where 


(85b) 
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PR*  =  Prandtl  number 


PR  =  C  y*/K* 
pa 


U*  =  viscosity  of  air  in  lb/ft-sec  evaluated  at  reference 
temperature 

K*  =  thermal  conductivity  of  air  in  Btu/ft-sec-°R  evaluated 


Therefore,  the  unit  surface  convective  heat  rate  for  high-speed  flow  (Ref.  16) 


VA  '  he  (TAK  -  T«) 


where 

qc/A  =  heating  rate  in  Btu/sec-ft* 

h  =  convective  heat  transfer  coefficient  in  Btu/ft* -sec-°R 
e 

A  =  local  wetted  area  in  ft-’ 

(1)  Heat  Transfer  Coefficient 

For  determining  the  heat  transfer  coefficient  in  laminar  flow  over 
a  flat  plate,  an  empirical  method  using  Blasius  theory  was  developed  bv  P.ubesin 
and  Johnson  (Ref.  17)  to  account  for  the  combined  effects  of  the  local  Mach 
number  and  the  temperature  ratio  of  the  wall,  T,^,,  to  the  boundary  layer  edge, 
TEMPE.  In  this  method,  the  conventional  analysis  of  heat  transfer  by  convection 
can  be  used,  but  the  reference  temperature,  T^^.,  for  evaluating  the  fluid 
properties  is  expressed  as  a  function  of  the  local  Mach  number  and  of  the  ratio 
T.,/TEM?E.  The  method  has  been  extended  to  cover  turbulent  flow,  and  has  proved 
to  be  both  convenient  and  useful. 

Therefore,  for  laminar  flow  the  reference  temperature  is 


tref  =  TE>5PE 


E  [1  +  0-58(fS?f-  l)  +  0-032  MACHt' 


and  for  turbulent  flow 


TEMPE  1  +  0.45  ^-,"pE  ~  l)  J-  0.033  MAC-HF. 


I 
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Equations  (87)  and  (88)  can  be  expressed  in  the  form  for  laminar  flow  as 

T^j.  =  TEMPE  +0.58  -  TEMPE^  +0.19  (TAW  -  TEMPE)  (89) 

and  for  turbulent  flow  as 

=  TEMPE  +  0.45  ^Ty  -  TEMPE j  +  0.21  ^TAy  -  TEMPER  (90) 

The  heat  transfer  coefficient  (Refs.  16,  18)  is  related  to  the 
flow  properties  through  the  expression 


1/3  / 

Nu* (PR*)  '  =  CfR  V 


(91) 


where 


N*  =  Susselt  number  =  h  x/K* 
u  e 

h  =  effective  heat  transfer  coefficient 
e 

X  =  distance  from  nose,  ft 

K*  *  thermal  conductivity  of  air,  Btu/ft-sec-°R,  evaluated  at 
T 

aR£F 

Rg  =  B/L  edge  Reynolds  number  =  P0v£X/ue 

oe  =  boundary  layer  edge  density,  lb-sec2/ft4 
=  boundary  layer  edge  velocity,  ft/sec 
C  =  0.575  for  laminar  flow  on  cone 
C  =  0.0296  for  turbulent  fiow  on  cone 
C  =  0. 778  for  laminar  flow  on  blunt  body 
C  =  0.0348  for  turbulent  flow  on  blunt  body 
a  *  0.8  for  turbulent  flow 
a  =  0.5  for  laminar  flow 

u  =  3/L  edge  viscosity,  lb-sec/ft‘ 


The  solutions  to  these  equations  for  the  neat  transfer  coefficient  is  obtained 
by  the  use  of  the  Blasius  incompressible  flat  plate  skin  friction  coefficients 
modified  for  compressible  flow  by  use  of  Eckert's  reference  enthalpy.  For 


48 


AFWL-TR-66-61 


t 

i 


laminar  flow 


and  for  turbulent  flow 


(92a) 


(92b) 


Good  correlation  with  experimental  data  is  obtain  .-d  if  the  gas 
properties  are  evaluated  at  the  reference  temperature  and  the  velocity  is  taken 
as  the  boundary  layer  edge  velocity.  For  the  case  of  mass  addition,  the 
properties  of  the  injected  gas  are  used  to  compute  the  Prandtl  number  and 
Nusselt  number. 


It  is  important  in  the  selection  of  a  suitable  material  for  thermal 
protection  by  the  ablative  process  that  the  produces  of  decomposition  have  a 
high  specific  heat.  This  in  turn  produces  a  high  effective  mean  specific  heat 
of  the  gas-air  mixture  in  the  boundary  layer  and  a  high  Prandtl  number.  It  is 
also  desirable  that  the  ablator  have  a  low  thermal  conductivity.  This  will 
decrease  heat  conducted  to  the  interior  of  the  vehicle.  The  surface  of  the 
vehicle  also  receives  heat  by  radiation  from  the  hot  gas  in  the  shock  layer  in 
addition  to  that  by  convection.  As  the  air  density  is  low,  its  emissivity  is 
much  less  than  1  percent.  Therefore,  the  rate  of  radiation  heat  transfer  is 
negligible  compared  with  that  of  aerodynamic  convection  and  can  be  ignored. 

b.  Unsteady  Heat  Conduction 

A  means  of  applying  the  Schmidt  graphical  method  for  solving  an  unsteady 
heat  conduction  problem  is  presented  here.  A  thorough  discussion  of  this 
method  can  be  found  in  any  general  text  (Refs.  16,  19)  and  will  not  be  attempted 
here.  This  method  is  quite  flexible  in  that  difficult  boundary  conditions  can 
be  handled  easily. 


This  method  is  widely  used  for  solving  unsteady  heat  conduct 
because  it  gives  an  iterative  profile  of  the  temperature  change.  5,. 
its  simplicity,  mistakes,  if  the>  occur,  are  easily  found,  ar.d  relat 
untrained  personnel  can  accomplish  the  work.  Numerical  methods  art- 
convenient  when  a  high-speed  digital  computer  is  available  since  the 
a  numerical  solution  car.  be  programmed  relatively  casiiv. 


ion  problems 
cause  of 
i  vt  i  v 

t special  ly 
steps  in 
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The  numerical  method  for  solving  unsteady-state  heat  conduction  differs 
from  that  used  to  solve  steady  state.  In  the  latter  case,  the  temperature 
distribution  in  a  body  can  be  obtained  for  a  network  of  points  in  a  solid  by 
solving  a  system  of  residual  equations.  In  unsteady-state  systems,  the  initial 
temperature  profile  is  known,  but  its  variation  with  time  must  be  determined. 
Therefore,  it  is  necessary  to  resolve  the  temperature  profile  at  some  future 
time  from  a  given  distribution  at  an  earlier  time. 

To  illustrate  the  numerical  method,  it  is  necessary  to  transform  the 
Fourier  conduction  equation  (a  partial  differential  equation  that  is  second 
order  in  space  and  first  order  in  time)  for  the  unsteady  temeprature  distribu¬ 
tion  in  a  heat-conducting  solid  into  a  finite  difference  form. 

A  T  aA2T 

_L_  =  -Ji_  (93) 

(Ax)2 

The  subscripts  denote  the  differentiation  variable.  Letting  n  denote  position 
and  t  time,  A_T  can  be  written  as 


A  T  =  T  -  T 

-  n,t+l  n,t 


and  in  a  similar  manner 


A  T  -  T  -  T  _ 

x  n,t+l  n,t 


The  exoression  A2T  thus  becomes 
x 


A2T  =  T  ,  -  2T  +  T  , 

x  n+l,t  *.i,t  n-l,t 


Substituting  these  expressions  into  the  Fourier  equation  (93),  gives 


_  _  _  aAt  t-r  -  2T  t  \ 

i,t-rl  *n,t  "  r .  \2  \  n+l,t  n,t  '  Si-1, tj 


The  temperature  throughout  a  wall  or  slab  can  nov  be  computed  for  any 
later  time  if  the  initial  distribution  is  known.  A  Schmidt  plot  demonstrates 
the  temperature  profile  versus  time  in  a  semi-infinite  slab  (figure  13). 
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Figure  13.  Schmidt  Plot  in  an  Infinitely  tnick  Wall 

A  constant  v«,ll  temperature  is  assumed  in  the  graphical  example,  but  a 
varying  wall  temperature  can  be  handled  with  equal  ease.  By  letting  the  wall 
temperature  vary  with  time,  the  distribution  at  each  point  within  the  body  can 
be  computed  for  each  time  increment. 

c.  Vehicle  Thermal  Model 

To  determine  the  time-temperature  profile  in  the  vehicle  skin.  tht. 
Fourier  conduction  equation  was  written  as  a  finite  difference  equation  and 
solved  numerically -  This  method  is  a  further  adaptation  cf  the  Schmidt  graph¬ 
ical  method  (Ref.  19).  The  wall  is  divided  into  a  number  of  lamina  with  known 
thickness,  specific  heat,  and  thermal  conductivity  and  with  a  known  initial 
temperature  distribution.  Figure  14  illustrates  the  model  with  its  a>sociated 
nomenclature. 
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Figure  14.  Model  for  Wall  Temperature  Profile  and  Ablation 
Recession  Calculations 


The  difference  equation  for  computing  the  temperature  at  position  n  at 
time  t  +  1  is 


T  =  T  + 

n,t+l  n,t 


aAt 

(Ax)2 


(Tirfl,t  ”  2Tn,t  +  Tn-l,t) 


(95) 


In  this  form,  the  material  properties  cannot  vary  from  layer  to  layer.  If  the 
properties  vary  from  lamina  to  lamina,  the  above  equation  must  be  written  in 
the  form 


AFWL-TR-68-61 


Where  a  is  the  thermal  diffusivity  and  may  be  replaced  by  its  defined  equivalent. 


where 


K  *  thermal  conductivity  of  material  in  Btu/ft-sec-°R 
p  *  density  of  material  in  lbs/ft3 
cp  *  specific  heat  of  material  in  Btu/lb-°R 


This  step  allows  one  to  easily  recognize  certain  groups  of  terms  as  the  heat 
conducted  through  each  lamina  and  permits  the  convective  heat  flux  to  be  used 
in  solving  for  the  wall  temperature.  The  equation  for  the  wall  temperature  is 


T  =  T  + 

n-l,t+l  n-l,t 


- -  [h  (Tx\W  -  T  \1 


2At  K 


(p  fix2) 

'  cp  rn-1 


(Tn-l,t  ~  Tn,t) 


The  backface  boundary  condition  can  be  specified  in  several  ways.  If 
internal  heating  <or  cooling)  is  present,  the  backface  may  be  held  constant  or 
allowed  to  vary  in  a  specified  manner.  A  conservative  assumption  is  that  no 
heat  flows  through  the  last  lamina.  This  assumption  will  ca-’se  somewhat  more 
ablation  and  higher  temperatures  than  would  be  encountered  with  a  cooled  back- 
face  or  other  heat  sink  material. 

d.  Ablation 

Advances  in  hypersonic  atmospheric  flight  have  resulted  in  environments 
of  extremely  high  temperatures.  Boundary  layer  temperatures  in  excess  of 
10,000°F  are  characteristic  of  vehicles  entering  the  atmosphere  at  hypersonic 
speeds. 

In  hypersonic  flight,  it  becomes  obvious  that  capacitance  or  mass  heat¬ 
sink  protection,  though  simple,  is  an  inefficient  means  of  contending  with  the 
extremely  high  heat  fluxes  associated  with  certain  reentries  or  with  sustained 
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periods  of  thermal  flight.  Consequently,  other  means  of  cooling  or  protecting 
for  operating  beyond  heat-sink  limitations  must  be  used.  One  convenient  means 
of  environmental  protection  and  for  which  this  program  was  primarily  written  is 
the  ablation  cooling  method. 

The  ablative  process  is  illustrated  in  figure  15.  The  ablation  process 
works  in  the  following  manner:  (1)  the  material  or  ablator  acts  as  a  heat 
sink;  (2)  when  the  critical  or  melting  temperature  of  the  ablator  is  reached, 
a  thin  layer  of  the  material  at  the  surface  will  begin  successively  to  melt, 
vaporize,  depolycerize,  or  decompose  chemically;  and  (3)  as  the  material 
vaporizes,  the  gaseous  products  of  decomposition  enter  into  the  boundary  layer. 
Being  relatively  cool,  as  compared  to  the  boundary  layer  air.  the  injected  gas 
forms  a  thin,  but  effective  film  that  reduces  the  heat  transfer  to  the  vehicle 
skin.  It  is  important  in  the  selection  of  a  suitable  material  for  thermal 
protection  by  the  ablative  process  that  the  thermal  conductivity  of  the  material 
should  be  as  low  as  possible  so  as  to  confine  the  high-temperature  zone  at  the 
surface  to  the  thinnest  layer  possible.  Likewise,  to  reduce  the  rate  of  mass 
loss,  and  hence,  the  required  weight  of  protective  ablator,  the  surface  tempera¬ 
ture  at  which  decomposition  and  melting  begins  should  be  as  high  as  practicable. 
The  products  of  decomposition  should  preferably  have  a  high  Prandtl  number, 
since  it  is  defined  as  ratio  of  heat  storage  to  heat  conduction  of  a  gas.  This 
means  it  is  desirable  to  have  the  mean  effective  specific  heat  of  the  gas-air 
mixture  in  the  boundary  layer  as  large  as  possible. 
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or  by  substituting  approximate  equations  for  each  gives 

1/- 


Q* 


e»(TS-Tb)  +  0-7(-^)  \  <“>, 


i1  -  ’rad^o) 


(102) 


where 


A 


Cp  »  heat  capacity  of  material 

*  ablative  surface  temperature 
»  temperature  of  unheated  material 

air  *  solecular  weight  of  air 

*  molecular  weight  of  injected  gas 
ov  «  rate  of  mass  injection 

(ih)Q  =  ..-.thslpy  difference  across  boundary  layer  without  t’vnspiration 

If  the  radiation,  Qra^»  Is  insignificant,  equation  (102)  can  be  reduced  to  the 
form 


Q*  -  A  +  B  (Ah) 


(103) 


where  A  is  an  empirical  value  obtained  from  experimental  investigations  of 
different  ablators  and  B  is  dependent  on  whether  flow  is  laminar  or  turbulent. 

Die  surface  recession  rate  is  obtained  by 

♦ 

S  -  jfc  (104) 

s  *  recession  rate,  ft/sec 
o  *  ablation  material  density,  lt-/ft3 

The  enthalpy  difference  across  the  boundary  can  be  calculated  from  the  following 
equations: 


where 


(Ah) 


h  -  Cp  T 
r  a  V 


(105) 
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where 


h  -  Cp  •  TEMPE  +  RF  •  V 
r  a 


I'M 


The  following  are  analytical  expressions  used  in  this  program  for  deternining 
mass  loss  and  recession  rates  of  three  common  ablators. 

(a)  Phenolic  Refrasil  (Silica) 

The  expression  for  obtaining  the  effective  heat  of  ablation 
for  Phenolic  Refrasil  is 


Q*  -  4710  +  3  (Ah) 


(106) 


where 


3  *•  0.58  for  lan  in  a  r  flow 


3  *  0.32  for  turbulent  flow 


(b)  ATJ  Graphite 

The  effective  heat  of  ablation  for  ATJ  graphite  is 


Q*  -  2000  +  3  (Ah) 


007) 


where 


3*2  for  either  turbulent  or  laminar  flow 


(c)  Carbon  Phenolic 


The  expressions  for  carbon  phenolic  are  core  complicated 
than  those  used  for  the  other  materials  due  to  reaction  of  carbon  with  atmos¬ 
pheric  oxygen.  This  material  has  a  very  high  melting  temperature  and  a  low 
thermal  conductivity  which  gives  it  excellent  insulation  properties.  Also,  its 
products  of  decomposition  have  a  high  specific  heat  compared  to  that  of  air- 


For  laminar  flow,  the  diffusion  mass  loss  ra.e  is  given  fev 


(Refs.  20,  21) 


XI,  +  X2,  (h  -  Cp_,  Tvl 
L  l  v  :  roL  / 


(IDS) 
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(2)  Total  Weight  Loss  Rate  and  Total  Weight  Loss 

After  the  local  instantaneous  and  total  local  recession  rates  have 
been  computed,  the  local  weight  loss  rate  can  be  calculated,  if  needed,  by 
determining  the  voltxse  reduction  of  the  local  segaent  and  multiplying  by  the 
material  density.  But  the  primary  weight  loss  of  interest  is  the  total  for  the 
vehicle.  To  compute  the  total  vehicle  mass  loss  end  mass  loss  rate,  the 
following  procedure  is  used.  Using  the  local  weight  loss  rate,  the  instantaneous 
and  total  recessions,  and  the  initial  radii  and  slant  length  of  the  particular 
segment  of  interest  (see  figure  16),  the  mass  loss  rate  per  unit  area  is 
multiplied  by  the  appropriate  surface  area  to  give  a  local  integrated  mass  loss 
rate.  This  is  accomplished  by  the  following  equations: 

iT(LES)  *  ?I  (^(LES)  +  ^(LES)  ~  2^(LEH>)  ^ c  S°L(LEK)  (112) 


where 


• _ 

'‘T(LES)  *  instantaneous  mass  loss  rate  of  segment  (LE5),  lb/sei 
(LEX)  -  number  of  particular  segment  of  cone 
W.  *  larger  segment  radius,  ft 
RS  *  smaller  segment  radius,  ft 
C(lejj)  *  total  recession  of  segment,  ft 

S  *  instantaneous  recession  of  segment,  ft/sec 
c>c  «  surface  material  density,  lb/ft3 
SOL  *  slant  height  of  conical  segment,  ft 


K  +  VT 


(LES) 


(113,* 


will  give  total  vehicle  ?  ->ss  loss  rate  when  the  operations  in  equation  (113) 
are  accomplished  in  a  program  "Do  Loop."  Summing  these  local  rates  will  result 
in  the  total  vehicle  weight  loss. 

2.  Stagnation  Heating 

To  evaluate  stagnation  point  convective  heat  transfer  rates,  this 
program  employs  Lee’s  theory  (Ref.  22),  modified  by  use  of  Eckert’s  reference 
enthalpy  techniques  (Ref.  19),  in  Lee's  equation 
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where 


q 


CSTAG 


se 


(114) 


7  —1.10  -  1.20  at  high  temperatures 

h  -  Y2/2CJ 
se  » 

T  *  1.4  ratio  of  specific  heats  for  air 


CWS 

Vvs 


dl’/dx 


density  evaluated  at  stagnation  reference  conditions 

viscosity  evaluated  at  stagnation  reference  conditions 

velocity  ccepeoent  gradient  at  the  stagnation  point 
(iron  Fay  and  Riddell) 


Fg  ■  stagnation  pressure  (oblique  shock) 

(1)  Sose  Blunting 

In  the  determination  of  nose  blunting,  it  is  assumed  that  the 
initial  nose  shape  of  particular  interest  was  a  sphere-cone,  and  after  ablation 
the  final  shape  was  a  sphere-cone-  Experimental  evidence  has  shown  that  the 
final  eroded  shape  could  he  approximated  by  a  sphere-cone  having  a  snail 
deviation  from  a  sphere.  The  reference  sphere-cone  configuration  after  erosion 
is  obtained  by  placing  a  spherical  surface  tangent  tc  a  cone  parallel  to  the 
original  surface  but  displaced  by  the  erosion  on  the  cone  and  is  located  axially 
at  a  position  on  the  vehicle  center  line  in  accordance  with  the  erosion  at  the 
stagnation  point  (figure  17). 

Tbe  experimental  evidence  of  nose-shape  change  of  reentry  vehicles, 
as  previously  noted,  can  be  represented  by  a  sphere-cone  which  deviates  from  a 
spherical  shape  by  an  asoun t  less  than  10  percent  of  the  final  nose  radius. 

The  final  eroded  sose  radios  can  be  obtained  in  terns  of  the  stagnation  point 
erosion  and  cone  erosion  as  shown  in  the  following  equation: 
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Figure  17.  Vehicle  Nose-3hape  Change 


The  integral  represents  the  erosion  at  a  given  time  during  the 
flight  for  the  locations  at  the  stagnation  point  S  and  a  point  on  the  cone  C  at 
a  location  of  wetted  length  of  about  five  times  the  nose  radius. 

It  is  now  possible  to  evaluate  the  final  nose  radius  for  a  given 
application.  Let 


(115) 
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and 


t 


(116) 


therefore 


R 


R  + 
o 


1 

(l  -  sin 


[sin  0 

L  c 


(117) 


where 

Rq  =  original  nose  radius,  ft 
R  =  instantaneous  nose  radius,  ft 

6c  *  cone  half-angle 

Ss  *  stagnation  point  recession  rate,  ft/sec 
C  «  recession  rate  at  location,  C,  ft/sec  (see  figure  17) 

f.  Boundary  Layer  Transition 

The  problem  of  transition  from  laminar  to  turbulent  boundary  layer 
fiow  has  always  been  a  troublesome  one  for  the  aeronautical  engineer.  The 
standard  approach  usually  has  been  to  design  conservatively,  that  is,  for 
turbulent  flow.  In  general,  the  transition  Reynolds  number  has  been  found  to 
primarily  depend  upon  the  local  Mach  number,  which  has  been  observed  in  the 
analysis  of  boundary  layer  transition  on  a  series  of  flight  vehicles.  Transi¬ 
tion  on  blunt  spherical  nose  vehicles  with  low  boundary  layer  edge  Mach  numbers 
appears  to  occur  at  Reynolds  number  of  1.5  x  106,  while  on  sharp  bodies  where 
the  edge  of  the  boundary  layer  edge  Mach  numbers  approach  10  Reynolds  numbers 
as  high  as  2.0  x  107  have  been  observed  prior  to  transition.  Thus,  for  untested 
vehicle  configurations,  it  is  necessary  to  investigate  the  relationship  of  the 
Mach  number  and  the  Reynolds  number  in  both  the  high  ar.d  low  Mach  number  regions. 
This  relationship  between  Mach  number  and  Reynolds  number  stands  to  reason  if 
their  definitions  are  understood.  The  parameter  Mach  number  describes  the 
influence  of  compressibility  on  heat  transfer  and  flow  phenomena  and  is  defined 
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as  the  ratio  of  the  gas  or  flight  velocity  to  the  local  or  ambient  speed  of 
sound.  The  Reynolds  number,  which  describes  the  nature  of  flow,  is  a 
dimensionless  measure  of  the  ratio  of  inertial  to  viscous  forces. 

It  should  be  noted  that  flight  data  on  sharp  bodies  indicate  that, 
generally,  transition  does  not  occur  instantaneously  over  the  entire  vehicle 
so  that  it  may  travel  several  thousand  faet  between  the  onset  of  transition 
and  tne  establishment  of  a  fully  turbulent  boundary  layer.  See  references  23 
and  24. 

The  approach  used  in  this  program  is  to  apply  the  sharp  and  blunt  body 
transition  Reynolds  numbers  simultaneously.  Therefore,  the  sharp  body  criterion 
is  described  by  the  expression 


ReTRAN  =  6'6  X  10? 

and  the  blunt  body  criterion  is  applied  when  either  of  the  following  two 
Reynolds  numbers  are  reached  with  conditions  stated 

ReTRAN  “  1,5  x  106  on  the  spherical  nose 

ReTRAN  *  5*°  x  106  on  the  conical  portion  aft  of  the  spherical 
nose  for  a  wetted  length  of  five  times  the  nose  radius. 

The  above  criteria  has  been  found  to  correlate  with  experimental  data. 
5.  AERODYNAMICS 

The  one  item  that  plays  the  most  significant  role  in  determining  the 
performance  of  a  high  velocity  or  hypersonic  reentry  vehicle  is  aerodynamic 
drag.  Experimental  results  have  been  correlated  with  theoretical  studies  of 
the  nonsteady  effects  on  a  reentry  vehicle  or  missile  during  its  trajectory. 
These  results  indicate  that  the  aerodynamic  forces  can  be  determined  accurately 
by  assuming  quasi-sbeady  flow  conditions;  that  is,  the  flow  field  at  any 
instant  of  time  is  the  same  as  that  associated  with  steady  motion  at  the  same 
velocity.  Furthermore,  it  has  been  demonstrated  that  the  boundary  layer 
behaves  in  a  quasi-steady  manner.  Therefore,  assuming  that  the  validity  of 
these  results  carry  over  for  an  ablating  blunt  cone,  the  instantaneous  drag 
coefficients  can  be  obtained  from  equivalent  steady-state  conditions. 
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This  program  asumes  that  the  drag  consists  of  three  components:  cone 
pressure  or  forebody  drag,  base  pressure  drag,  and  skin  friction  or  viscous 
drag.  It  is  convenient  to  consider  these  quantities  as  distinct  quantities 
which  can  be  added  together  to  obtain  the  total  vehicle  drag.  Though  all 
these  quantities  are  distinct,  some  are  dependent  upon  the  boundary  layer 
condition.  For  example,  the  condition  of  the  boundary  layer,  laminar  or 
turbulent,  which  influences  the  viscous  drag  can  also  influence  the  base  drag. 

a.  Cone  Pressure  Drag 

The  forefcody  pressure  drag  for  a  slender  body  at, zero  angle-of-attack 
or  at  incidence  can  be  calculated  on  the  basis  of  the  modified  siender-body 
theory. 

The  forebody  pressure  drag  is  defined  by  the  relation 


D  =  - 
P 


P  cos 


d  S 


m 


where  cos  ^n,  V^j  is  the  cosine  of  the  angle  between  and  the  outward  normal 
to  the  vehicle  sui 
except  base  area. 


to  the  vehicle  surface.  The  area  comprises  the  total  area  of  the  vehicle 


For  a  sharp  cone  at  a  *  0,  the  pressure  P£  is  constant  over  the  entire 
body,  and  the  forebody  drag  is  then  only  a  function  of  the  pressure  ratio 
P  /P^  across  the  conical  shock.  Solutions  for  this  pressure  ratio  as  a  function 
of  Mach  number  and  cone  angle,  8  ,  have  been  tabulated  by  Kopal  (Ref.  25)  for 
the  ratio  of  specific  heats  of  y  =  1.405.  These  r  suits  have  been  employed 
directly  in  the  drag  equation  in  order  to  obtain  the  forebody  drag  on  a  sharp 
cone.  The  analysis  produced  the  equation 


4  sin2  3 

c 


2.5  +  8  M  sin  9 
_ * _ c 

1  +  16  M  sin  r 
*  c 


(118) 


which  is  similar  to  tha  form  of  the  Newtonian  drag  equation  for  a  sharp  cone. 
The  last  term  represents  the  variation  from  the  Newtonian  theory  due  to  Mach 
number  effects. 
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Nose  blunting  of  a  cone  will  increase  the  pressure  drag.  Local  over¬ 
pressures  are  induced  near  a  blunt  nose  by  the  strong  nose-shock  curvature. 

The  results  of  sharp  and  blunt  body  solutions  have  been  correlated  as  a  function 
of  6c,  and  R^/Rg.  However,  for  bluntness  ratios,  Rj./Rg  less  than  0.05, 
the  values  for  the  drag  due  to  the  overpressures  are  low  compared  to  . 

This  program  do.  ^  account  for  the  drag  due  to  the  nose  being  blunt 
within  the  limits  expressed  above.  The  sharp  cone  drag  is  modified  by  assuming 
the  nose  drag  to  be  that  of  a  sphere  minus  that  of  the  cone  replaced  oy  the 
sphere.  Because  both  of  these  terms  are  referenced  to  the  projected  area  of 
the  sphere,  the  area  or  bluntness  ratio,  Rjj/Rg,  term  is  needed  as  a  correctness 
factor.  Therfore,  the  pressure  coefficient  for  a  blunt  cone  is  (Ref.  26) 


Dp 


cone 


-Dp 


(119) 


where  0.58  is  drag  coefficient  for  a  sphere  obtained  from  reference  27 
CDp  is  from  equation  (118) 

Rjj/Rg  is  the  bluntness  ratio  of  nose  radius  to  cone  base  radius 
b.  Skin-Friction  Drag 

The  skin-friction  or  viscous  drag  is  defined  as 


where  t  is  the  local  skin  friction  per  unit  area  due  to  viscosity,  and 
cos  ^t,  is  the  cosine  of  the  angle  between  and  the  tangent  to  the 
vehicle  surface  in  the  t  direction  as  shown  in  figure  18. 

The  skin-friction  drag  coefficient  can  be  defined  as 


C 


Df 


L 

* 

cos  r  dx 

ct 
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Figure  18.  Aerodynamic  Body  Subject  to  Normal 
and  Tangential  Forces 


where  is  the  local  skin-friction  coefficient  defined  as 

CO 


C  =  - - -  (121) 

-  1/2  o  Vj 

(1)  Local  Skin-Friction  Coefficient 
(a)  Laminar  Flow,  a  =  0 

The  local  skin-friction  coefficient  is  calculated  using  the 
Blasius  flat  plate  in  compressible  solution.  The  flat  plate  solution  was 
modified  for  conical  flow  by  the  Mangier  transformation  (Ref.  28)  and  for 
compressibility  by  Eckert's  reference  enthalpy  method  (Ref.  29). 

The  modified  Blasius  equation  is 


where 


Re  *  — - -  boundary  layei  edge  Reynolds  number 

X  u 

e 

x  *  distance  along  sharp  cone  surface 


022) 
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og  =  density  of  air 
ue  -  viscosity  of  air 
*  velocity 

The  superscript  *  and  subscripts  e  and  »  signify  the  property  is  evaluated  at 
the  reference  enthalpy  or  conditions  at  the  edge  of  boundary  layer  or  free 
stream,  respectively. 

Equation  (121)  can  be  simplified  to 


1.15 


V-1.5 

e 


P 


SC  CD 


(123) 


Letting 


P  V 
OD  00  X 


and  rearranging  equation  (123), 


1.5 


From  references  30  and  31,  the  quantity 


can  be  replaced  by 


where 


(124) 
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It  is  assused  that  cone  pressure  pc  is  equal  to  boundary  layer  edge  pressure, 
Pe»  the  free-strean  specific  heat  is  equal  to  the  boundary  layer  edge  specific 
heat,  the  specific  heat  ratio  is  y  *  1.405,  and  the  recovery  factor,  R^,  is 
equal  to  0.8426.  Therefore,  equation  (123)  for  the  local  skin-friction 
coefficient  can  be  written  as 


,  1.5  0.5  ,  .-0.185 

ipft)  ft)  ft)  ft)' 


-0.185 


(125) 


where 


h* 

—  -  0.5  +  0.5 

n 

e 


eft)’  *  • 


0374  Mz 
e 


(126) 


The  ratios. 


ft) 


ft) 

ft) 


are  obtained  frora  conical  flow  results  given  by  Bertran  (Refs.  31  and  32) 
which  have  been  correlated  as  a  function  of  the  hypersonic  sinilarity  parameter, 
sin  Sc  *  Kc»  A  curve-fit  of  these  results  produced  the  following  relations: 


t-  1  -  tr  W 


(127) 


—  *  1  +  2.8  K‘ 
?  c 


2.5  +  8  K 
_ c 

1  +  16  K 


(128) 
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T 

~  *  1  +  0.0966  Kc  +  0.2267 


(129) 


The  local  Mach  number  can  be  computed  from 


Y_R  r 


•fST 
e  e  e 


Froa  previous  assumptions  and  reducing 


*-  v-  \rJ 


(130) 


The  local  laminar  skin-friction  is  now  defined  in  teres  of  f re e-s cream  conditions, 
vail  temperature,  and  cone  angle. 

(b)  Turbulent  Flow,  a  *  0 

Using  the  Blasius  flat  plate  incompressible  solution  modified 
in  a  similar  manner  lor  conical  flow  and  coepressibility,  the  Blasius  equation 
for  the  local  skin  friction  coefficient  in  turbulent  flow  is 


2  D  V2 
e  e 


(131) 


By  substituting  in  a  manner  similar  to  laminar  flow,  and  using  the  same  ratios 
and  rearranging,  the  following  expression  for  the  turbulent  skin-friction 
coefficient  in  terms  of  free-stress  conditions,  vail  temperature,  and  cone 
angle  is  given  by 


(132) 


where 


0.038S  M- 
e 


for  Rp  =  0.877. 
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(2)  Average  Skin-Friction  Coefficient 

To  obtain  the  average  skin-friction  coefficient  for  both  laminar 
and  turbulent  flow,  the  local  value  must  be  integrated  over  the  cone  surface. 
Therefore,  substituting  -R2  for  in  equation  (120)  and  rearranging 


2- 


* 

a 

I 


C^.  cos  6  (r)  dx 


(133) 


but 


and 


-kur 


R2  -  L2  sin2  ec 


r  *  x  sin  S 


Substituting  these  relations  into  equation  (133),  then 


Df  L2  sin2 


-  \  (L)  cos  5^  sin  3^  |  (x) 

“c  x=L  -L 


dx 


(134) 


Kew,  combining  equation  (134)  with  equations  (125)  and  (133),  the  average 
skin-friction  coefficients  for  laminar  and  turbulent  flow  are,  respectively. 


ft)  'ft)  er 


and 


ftj  ft)  (W‘ 


(135) 
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The  wall  teaperature  used  in  the  enthalpy  ratio  equation  is  taken  as  an 
average  teaperature  over  the  entire  vehicle  surface. 

c.  Base  Drag 

The  base  pressure  drag  is  determined  by  the  mechanics  of  the  wake,  for 
which  there  is  as  yet  no  complete  theory.  Values  of  the  base  pressure 
coefficient  Cpg  must  be  obtained  experimentally.  Because  of  this,  analytical 
methods  of  treating  base  drag  are  usually  replaced  by  a  correlation  of  experi¬ 
mental  results.  As  the  vehicle  velocity  increases,  the  base  pressure,  and 
therefore,  the  base-drag  coefficient  are  affected  similarly  to  other  pressures. 
Once  the  speed  is  well  into  the  supersonic  region,  most  base  pressures  for  a 
cone  appmach  about  70  percent  cf  a  vacuus.  Hoerner  (Ref.  33)  has  correlated 
a  wide  variety  of  data  for  conical  flow  with  the  following  result: 

CpB  -  1.43/M*  (136) 

for  a  vacuus.  Reducing  this  by  70  percent 

CpB  -  1.001/U^  (137) 

for  M  _>  2. 

d.  Stability  Derivatives 

The  stability  derivatives  used  in  this  prograe  are  for  a  sharp  cone  and 
are  easily  derived  from  references  34  and  35.  These  derivatives  consisted  of 
the  damping  coefficient  per  radian  of  pitch  angle  cf  attack,  CKQ,  and  the  rate 
of  change  of  the  normal  force  coeficient  per  degree  of  pitch  angle  cf  attack, 

CS  .  These  values  are  for  the  following  cone  conditions: 

4.o*  <  e  <io* 
c 

and 

0  <  Rjj/Rg  <  0.3 

The  rate  of  change  of  pitching  soncnt  coefficient,  CKA,  can  be  computed 
as 

OJA  *  —  (SK)  (CSA)/  DIA  (OS) 
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where 

SH  *  static  aargin  in  feet 
OIA  -  reference  diameter  in  feet 

Because  of  the  symmetry  of  the  conical  vehicle  used  in  this  program,  the  values 
produced  by  yaw  angle  of  attack  and  yaw  angular  velocities  are  equal  in  magni¬ 
tude  to  pitch  (data).  Therefore, 


CSR  - 

CKQ 

(139) 

CSS  - 

-OiA 

(1*0) 

CT3  - 

-CSA 

(1*1) 

The  derivatives  are  used  to  compute  the  aerodynamic  forces  and  moments  needed 
tc  restore  the  vehicle  to  aero  angle  of  attack  in  pitch  and  yaw. 
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SUBROUTINES 


1.  ATMOSPHERIC  PROPERTIES 

The  following  list  of  subroutines  are  incorporated  in  STRAB-6  to  compute 
the  necessary  atmospheric  properties.  This  atmosphere  uses  the  equations 
derived  for  the  1962  U.S.  Standard  atmosphere. 

a.  Subroutine  RHOF 

This  subroutine  computes  the  atmospheric  density  as  a  function  of 
geometric  altitude  from  the  main  program.  The  geometric  altitude  is  converted 
to  geopotential  altitude  in  the  subroutine  and  the  density  is  computed 
accordingly.  The  units  of  density,  Rho,  are  lbs-sec2/ft4. 

b.  Subroutine  VSD 

This  subroutine  calculates  the  velocity  of  sound,  VS,  in  feet/sec  as 
a  function  of  geometric  altitude. 

c.  Subroutine  TEMPA 

The  ambient  temperature,  TEMP,  is  calculated  in  this  subroutine  as  a 
function  of  geometric  altitude.  The  units  are  degrees  R. 

d.  Subroutine  VISCO 

This  subroutine  computes  the  viscosity  of  air  as  a  function  of 
geometric  altitude.  The  units  are  lb-sec/ft2. 

e.  Subroutine  PRESS 

The  ambient  pressure,  PE,  is  calculated  in  this  subroutine  as  a 
function  of  geometric  a.  itude.  The  units  are  lb/ft2. 

2.  HOT  AIR  PROPERTIES 

a.  Subroutine  VISCOT 

This  subroutine  calculates  the  viscosity  of  air  as  a  function  of 
temperature.  It  is  used  primarily  for  the  hot  air  in  the  nonablating  boundary 
layer.  The  units  are  lb/sec-ft. 
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b.  Subroutine  TCAT 

The  thermal  conductivity,  TCA,  of  hot  air  is  calculated  ir  this 
subroutine  as  a  function  of  temperature.  The  units  of  TCA  are  Btu/ft-sec-°R. 

3.  ABLATOR  THERMAL  PROPERTIES 

The  following  subroutines  are  used  to  compute  the  thermal  properties  of 
carbon  phenolic,  phenolic  silica  (Refrasil),  and  ATJ  graphite.  The  properties 
computed  are  the  thermal  conductivity  and  the  specific  heat.  The  carbon 
phenolic  and  phenolic  silica  ablator  have  both  virgin  and  char  properties  as 
well  as  densities.  From  all  available  data  for  ATJ  graphite,  it  appears  that 
no  distinction  is  made  between  virgin  or  char  conditions.  The  subroutines  that 
compute  the  thermal  properties  for  carbon  phenolic  and  phenolic  silica  have  the 
same  names  to  minimize  card  changes  in  the  main  program. 

a.  Subroutine  VTCPC 

This  subroutine  computes  the  virgin  thermal  conductivity  for  carbon 
phenolic  or  for  phenolic  refrasil  as  a  function  of  wall  temperature.  The 
units  are  Btu/ft-sec-°R. 

b.  Subroutine  VCPPC 

The  virgin  specific  heat  for  carbon  phenolic  or  phenolic  silica  is 
calculated  by  this  subroutine  as  a  function  of  wall  temperature.  The  units 
are  BtU/  lb-°R. 

c.  Subroutine  CTCPC 

The  char  thermal  conductivity  for  carbon  phenolic  or  phenolic  silica 
is  computed  by  this  subroutine  as  a  function  of  wall  temperature.  The  units 
are  in  Btu/ft-sec-°R. 

d.  Subroutine  CCPP 

This  subroutine  computes  the  char  specific  heat  of  carbon  phenolic  or 
phenolic  silica  as  a  function  of  wall  temperature.  The  units  are  Btu/lb-°R. 

4.  ABLATION  GAS  PROPERTIES 

Subroutines  GCPB  and  GASTC  compute  the  thermal  properties  of  the  ablative 
gas-air  mixture  in  the  boundary  layer.  These  subroutines  are  used  with  the 
carbon  phenolic  ablator  only.  They  compute  the  gas  specific  heat  and  gas 
thermal  conductivity.  The  values  associated  with  phenolic  refrasil  and  ATJ 
graphite  are  inputted  as  single  values. 
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5.  STABILITY  DERIVATIVES  SUBROUTINES 


These  subroutines  are  incorporated  into  STRAB-6  for  computation  of  the 
necessary  aerodynamics  discussed  in  section  II-5.  The  equations  are  general 
for  the  following  cone  parameters: 

4.0°  <  9  <  10.0° 

c 

0  <  V»3  <  0-3 
3.0  <  M  <30 

CO 

a.  Subroutine  CDSF 

This  subroutine  calculates  the  aerodynamic  drag  coefficients  CP,  CAN, 
CPB,  and  CDF  as  a  function  of  M  or  Reynolds  number.  CDSF  is  called  for  in  the 
main  program. 

(1)  Subroutine  CDEM 

This  subroutine  computes  the  skin-friction  drag,  CDF,  and  is 
called  into  the  program  through  subroutine  CDSF.  The  boundary  layer  edge 
properties  are  calculated  in  this  subroutine. 

b.  Subroutine  ACMQ 

The  damping  coefficient,  CMQ,  is  calculated  in  this  subroutine  as  a 
function  of  free-stream  Mach  number.  CMQ  is  per  radian. 

c.  Subroutine  ACNA 

This  subroutine  computes  the  normal  force  coefficient,  CNA,  as  a 
function  of  free-stream  Mach  mmber.  CNA  is  per  radian. 

6.  SUBROUTINE  INITIAL  (Ref.  36) 

This  subroutine  is  incorporated  into  STRAB-6  to  provide  the  necessary 
vehicle  geometry  calculations  and  changes  in  units.  By  reading  in  vehicle 
axial  stations  of  interest,  the  subroutine  will  compute  the  transition  Reynolds 
number,  slant  height  i.SOL) ,  smaller  radius  (RS),  and  larger  radius  (RL)  for 
each  conical  segment.  It  will  also  compute  the  appropriate  cone  base  radius 
and  total  slant  height.  See  figure  16  for  further  clarification. 


7.  SUBROUTINE  ATTK  AND  SSLP 


These  subroutines  are  used  to  compute  the  complete  range  of  angles  of 
attack,  ALP,  and  sideslip,  BFT.  These  angles  are  computed  as  a  function  of 
the  components  of  the  relative  velocity. 

8.  SUBROUTINE  DLON 

This  subroutine  computes  the  longitude  of  the  vehicle  with  respect  to  the 
earth  reference  system.  See  section  II-2c  for  derivation. 
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SECTION  IV 
PROGRAM  OPERATION 

1.  INPUT  DATA  REQUIREMENT 

All  input  data  that  is  subject  to  change  either  dje  to  vehicle  or  trajectory 
requirements  are  read  into  STRAB-6  from  the  data  statement.  The  word  data  is 
punched  first,  then  the  FORTRAN  symbol,  and  then  the  numerical  values.  The 
symbols  and  numerical  values  are  punched  in  Columns  7  through  72,  inclusive, 
with  up  to  nine  continuation  cards  if  needed.  Columns  73  to  80,  inclusive, 
can  be  used  for  comments  or  identification.  The  numerical  values  can  be 
punched  in  any  format  except  "I."  Care  should  be  taken  to  keep  the  symbol  and 
its  corresponding  value  in  the  same  sequence.  For  example: 

DATA  V.  AZN,  NOFX,  NAL,  GAMMA/22414.,  123.003,  8,  4,  -39.92/ 

Also, 

DATA  WE,  REO,  ET/0.729E-04,  20,927491E+C6,  8.15133302E-O2/ 

Another  method  of  data  card  is 

DATA  (POS  (I),  1=1,  6)/3. 338,  15.189,  27.040,  38.851,  50.742,  65.556/ 

Another  flexibility  of  the  input  routine,  aside  from  the  fact  that  the 
data  can  be  introduced  in  any  order  of  the  symbols  (with  their  associated 
values),  is  the  feature  that  permits  inputting  the  same  symbol  (and  an  associ¬ 
ated  value)  again,  thus  superseding  the  previous  value.  This  feature  is 
convenient  when  basically  constant  vehicle  configurations  and  initial  conditions 
are  maintained  in  many  runs  with  only  a  few  quantities  changing.  The  "Standard" 
input  data  cards  may  be  kept  intact,  and  only  the  varying  quantities  nay  be 
punched  on  a  card  after  the  data  cards  which  may  be  removed  later. 

For  the  heating  portion  of  the  program,  there  are  several  "Do  Loops"  that 
initialize  temperatures,  ablator  lamina  thickness,  control  indicators,  and 
thermal  properties  of.  the  adhesive  bond  and  backface  aluminum  shell.  Also, 
there  are  various  equations  to  obtain  initial  trajectory  conditions  as  well  as 
initial  dccamy  program  starters. 
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a.  Multi-Case  Run 

This  program  readily  lends  itself  for  a  multi-case  run  with  a  minimal 
change.  If  the  need  arises  for  a  parametric  study  of  various  vehicle  geometry 
and  trajectory  missions,  the  following  changes  can  be  made  to  some  of  the  data 
inputs  and  check  for  completion  of  the  multi-case  run. 

Read  1500,  NCASE,  (POS(J),  J=l,  8)  1500  FORMAT  (12,  7A10,  A8) 

Print  1501,  NCASE,  (POS(J),  J-l,  8)  1501  FORMAT  (1H1,  30X,  7A10,  A8) 

These  read  and  print  statements  allow  the  identifier  to  be  read  into 
the  program  and  case  heading  to  be  printed  out  before  that  particular  case  data 
output.  The  1500  format  allows  the  case  number  to  be  punched  in  the  first  two 
columns.  When  a  zero  is  in  these  columns,  the  program  will  go  to  end  with  this 
check  statement. 

If  (NCASE. NE.0)  go  to  1700 
END 

where  statement  1700  would  be  the  return  at  the  beginning  of  the  program  to 
sta^t  another  case. 

To  update  the  vehicle  configuration  and  trajectory  parameters,  the 
following  data  input  is  used 

Read  1300,  (POS(I),  I=l,NOFX) 

Read  1300,  THETC,  GAMMA,  BI 

Read  1300,  V,  AXN,  GLAT,  DLON 

Read  1300,  ALGTK,  SNR  1300  Format  (8F10.2) 

All  other  data  input  cards  would  remain  the  same.  See  appendix  II  for  multi¬ 
case  listing  and  setup. 

2.  OUTPUT  DATA  REQUIREMENTS 

The  output  quantities  of  the  program  and  sample  case  are  shown  in  appendixes 
I  and  III. 

Before  the  printout  of  any  results,  listing  of  the  program  and  input  data 
cards  are  printed.  This  record  of  the  actual  program  listing  and  data  input 
will  assist  in  the  identification  of  a  run  as  well  as  an  aid  in  troubleshooting 
in  the  event  of  arithmetic  errors. 


79 


hi . . .  i  t' 


t 


AFWL-TR-68-61 

The  first  block  of  data  printout  of  the  matrix  of  temperatures,  lamina 
thickness,  etc.,  are  obtained  from  the  initialization  "Do  Loops"  for  heating 
portion  of  the  program.  This  printout  is  done  only  once. 

The  second  block  of  data  is  the  trajectory  output,  aerodynamic  coefficients, 
and  boundary  layer  edge  values.  The  next  block  of  output  is  the  wall  and  skin- 
profile  temperatures,  total  recession,  integrated  weight  loss  rate,  and  total 
weight  loss. 

Printout  is  regulated  by  a  time  indicator.  Data  can  be  printed  out  at  any 
time  step  desired.  If  a  printout  of  data  is  wanted  every  0.1  second,  set 
BT  =  .1.  The  following  method  will  then  allow  printout  at  that  time  step. 

BTL  =  BT  -  H/6. 

BTH  =  BT  +  H/6. 

If  (T  .GT.BTL  .  AND.  T.  LT.  BTH)  go  to  (statement  ntanber) 

BT  =  BT  +  .1 

where  H  is  the  Runge-Kutta  integration  time  step.  This  method  will  allow  data 
printout  and  also  go  into  the  heating  portion  of  the  program. 

In  addition  to  the  printed  output,  selected  quantities  may  also  be  obtained 
as  punched  output.  The  punched  cards  may  be  useful  for  machine  plotting  or  as 
input  to  some  other  program. 
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PROGRAM  5TRAB62  < INPUT. OUTPUT) 

C  6  0E6REE  TRAJECTORY  WITH  ABLATION  HEATING  AND  STAG  A  2 

COMMON  /A/  HE  A3 

COMMON  /8/  AREA. RAD. SLANT. SNR. THETC.V  A  4 

COMMON  /C/  CPM.CP.CAN.CPB.COE  A  S 

COMMON  /O/  Tw(5o.12).IND{50.12).JIND(12).VCP(50.12).vTC{50.12).VRH  A  6 

10 (5n. 12), OEL (50. 12) »DE (50.12) »SD0T(12) .XI (12) .RS(12) ,RL(12) .SOL (12  A  7 

2)  ,WT(12) *TM { 12) «TWT(12) .SACP(SO) *SATC(50) .SARHO(SO) «RN(12>  *0STR(12  A  A 

3)  .ODOTC { 12) »QTOT ( 12) .QCONO ( 12)  A  9 

common  /£/  nofxm.nal.alt.amach.od.amache.tempe.rhoe.ve.visg.pbe.an  A  10 

IRE  A  11 

COMMON  /F /  NOFX.ALGTH.RAOI.POS(ll) .azn.glat  A  12 

VEHICLE  OAT*  A  13 

IF  NOFX  IS  CHANGED  CHANGE  FORMAT  103  A  14 

IF  CONE  ANGLE  IS  CHANGEO  SUB  CDSF  MUST  BE  CHANGEO  A  15 

DATA  SNR, THETC.BI.SM/. 0208. . 07853,6. 27836. .59718/  A  16 

DATA  PI. E.AJ,flRAV,G/3. 14159. 2. 7183. 778. .32. 174, 32-174/  A  17 

DATA  REO. PE. ET.Wt/20.927491E*06..673852E-02». 08181333. 0.7292E-04/  A  18 

DATA  NOFX. NAL.NALC.NALS/8. 4. 10*35/  A  19 

DATA  CELT.OT .H.RT. AT/.05* .0005* .001 . .05. 1 ./  A  20 

DATA  AZN.GLAT.OLON, GAMMA/154. *33.996.-107.5.-22.39/  A  21 

OATA  V. 2, ALPHA/15167., 2. 50E-*03. 7./  A  22 

DATA  (POS(l) .I*1.8)/6. 607*18,570. 30. 533*42,496. 54. 459. 66. *22. 78. 38  A  23 

15.86.469/  A  24 

*12*0.51  A  25 

ATX*23.o8  A  26 

CALL  INITIAL  ( SNR. THETC. SLANT)  A  27 

AL6THI*aLGTH  A  28 

SNR1«SNR  A  29 

CLAT1*ATAN((1,-ET**2)*TAN(GLAT))  A  30 

CLAT«CLAT1*1A0./PI  A  31 

ETA1.GLAT-CLAT1  A  32 

PHN* (90. -GAMMA-ALPHA) 4PI/180,  A  33 

PPHaO.O  A  34 

PH*(1 ,*PPH)*PHN  A  35 

DLGGaDLON  A  36 

PHlaPH  A  37 

VMAcaBI  A  38 

RAObRAOI  A  39 

N0FXP*N0FX*1  A  *0 

N0FXM*N0FX-1  A  *1 

AKL1*S370,  A  42 

AKTlc*240.  A  *3 

AKL?*5.37  A  44 

AKT2*5.77  A  *5 

GV I S«* • *?9E-5  A  *6 

TOL*. 00854  A  *7 

TMELT»6760,  A  48 

TEM*544.  A  *9 

NALTvNAL  A  50 

NP1*NAL*1  A  51 

NP2*NAL*2  A  52 

NALC1*NALC*1  A  53 

NALC2»NALC*2  A  54 

NPIS»NALS*1  A  55 

NP2S*NALS*2  A  56 

RHCV*90.4  A  57 

RhCC*7* .  A  58 

ATJ«»H0*l2ft.96  A  59 

CPAa.2395  A  60 

nn  1  Jal.NOFxP  A  61 
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JJND(J>«0 

A 

62 

SOOT ( J) *0. 

A 

63 

00  2  J*1 .NOFXM 

A 

64 

VCPtNPl.J)*.315 

A 

65 

VCP(NP2.J)*.208 

A 

66 

VTC(NPl.J)*.74E-04 

A 

67 

VTC (NP2.J) *.03694 

A 

68 

VRH0{NP1,J)*91.7 

A 

69 

VRM0(NP2. J)«169. 

A 

70 

00  2  1*1 » NP2 

A 

71 

TW ( I . J) *T£M 

A 

72 

IN0(I« J)*0 

A 

73 

DO  3  1*1. NAL 

A 

74 

00  3  J*1 .NOFXM 

A 

75 

DEL ( I . J) *TOL 

A 

76 

OE  < I »  J) *T0L 

A 

77 

00  4  1*1. NOFXM 

A 

78 

OEL(NPl.I)*. 00333 

A 

79 

DE<NP1»I)«. 00333 

A 

80 

DEL(NP2.I)*.005 

A 

Ml 

•?E<NP2.I)*.005 

A 

82 

00  5  1*1 «KALS 

A 

83 

INO l I »NOFXP>  *0 

A 

84 

oelci.nofxpj.tol 

A 

85 

OE  1 1 .NOFXP)  *T0L 

A 

86 

TEMP-299. 2 

A 

87 

T4M.TEMP 

A 

88 

PHl-80.967 

A 

89 

STAM-TEMP 

A 

90 

stol-tol 

A 

91 

SX*.01745*RS (NOFXP) *PHI 

A 

92 

RX*SNR 

A 

93 

WCP (NAUCl .NOFX) «.3l5 

A 

94 

VCP (NALC2. NOFX) « .208 

A 

95 

VTC (N4LC1 .NOFX) ..74E-04 

A 

96 

VTC (NALC2.N0FX) a.03694 

A 

97 

VRHO (NALC1 .NOFX) *91 ,7 

A 

98 

VRH0(NALC2. NOFX) *169. 

A 

99 

00  6  1*1  .NALC2 

A 

100 

TN(I,N0FX)»TEM 

A 

101 

1ND«!»N0FX)*0 

A 

102 

nEL ( I .NOFX) *TOL 

A 

103 

OE ( I .NOFX) *TOL 

A 

104 

DEUNALC1  .NOFX)  *.00333 

A 

105 

OE (N*LC1 .NOFX) *,00333 

A 

106 

DEL (NALC2.N0FX) *.005 

A 

107 

OE  CNALC2.N0FX) *,005 

A 

108 

nE (NP2S.N0FXP) *,005 

A 

109 

DEL (NP25*N0FxP) *.005 

A 

110 

OE (NPlS.NOFXP) *.00333 

A 

111 

OEL<NP1S.NCFXP)..C0333 

A 

112 

SACP(NPlS) *.315 

A 

113 

SACP(NP2S)*.20A 

A 

11* 

SATC<NPlS)*.74E-04 

A 

115 

SATC (NP2S) *.03694 

A 

116 

SARHO (NP1S) *91.7 

A 

117 

SARHO  <NP2S) *169, 

A 

118 

00  7  N*1 ,NP2? 

A 

119 

Tu (N.NOFXP) *TEM 

A 

120 

PSI-O. 

A 

121 

0AZ*0.0 

A 

122 

AZ«(1.000*OAZ)*AZN 

A 

123 
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SGL*T«SIN<GL4T) 

CGLATaCOS(GLAT) 

SAZ»SIN(AZ) 

C4Z«C0S<*Z) 

W£X«»W£«cgLAT*SAZ 

w£r«WE»CGLAT«CAZ 

kEZ»wE»SGLaT 

Rei!£!Q/tl*!,00*PE#tSIN<CUT1>>##2,#**500 

B«H,P£1#SIN(£TA1J*SAZ 
rYH«-RE]«S!N(ETa1)*CaZ 
RZH»REl*Cns<ETA} ) 

T*0,0 

*■0.0 

DX*0*0 

Y«0.0 

PPHaO.O 

PM*(1 «0*PPH) *PHN 
0Z«V*SINt0*MMA*PI/l80.) 

DY»V*COS (GAMMA»Oi/ieo. ) 

LLL«0 

alt«z 

TH«0,0 

DTH«0,0 

OPHan.O 

PS«C.O 

OELPSaO.O 

DPSNaP.apf 

0P5« { 1 . ♦OELPS) »DPSN 

U«0.U07639E*17 

BJ»Ift23.*JE-fl6 

BH»6.o*E-06 

8Kafe.37E-06 

A9EA*PI*R*p**2 

OIA.i.781 

DCDaO.OO 

D2Xa0*0 

DPYaO.O 

02ZaO.O 

OPTHaO.O 

02PHaO.O 

02PS»0.0 


PRINT  71 
PRINT  72, 
PRINT  9?, 
PRINT  73, 
PRINT  7*, 
PRINT  75, 
PRINT  74, 
PRINT  77, 
PRINT  78, 
PRINT  79, 
PRINT  71 
PRINT  an 


<JINO(lJ ,Ial,NOFXP) 
fSDOT(I) ,Ial,NOPXP) 
l  lINDd,  J)  , Jal «NOFX)  ,Iai,NALC2) 
««VCP(I.JJ :J»l,NOFX},Ial,NALC2) 

{ <VTC(I, J) *J*l,NOFX) ,I«1,NALC?) 

{ (VRHO ( I »  J) , Jaj ,NOFX) , laj , NALC2) 
C COEL { I , J» , Jal ,N0FX) , lal ,NALC2 J 
t  <OE(I,J) , J*1 »NOFX) ,Ial,NALC2) 

( (TN { T , J) , Jal «NOFX) , lal «NALC2) 


CALL  V50  (VS, ALT) 
CALL  RHOF  (RHO.ALT) 
AMXCHaV/VS 


CALL  VISCO  (VIS, ALT) 
ANRaRHO*V*SL ANT/VIS 
CALL  COSE  (CO, AmACH, ANR) 
Qn«,5*R«0*V*»? 
ALPToABS(ALPMA) 


A  124 
A  125 
A  124 
A  127 
A  128 
A  129 
A  13o 
A  131 
A  132 
A  133 
A  134 
A  135 
A  136 
A  137 
A  138 
A  139 
A  140 
A  141 
A  142 
A  143 
A  144 
A  145 
A  146 
A  147 
A  148 
A  149 
A  150 
A  151 
A  152 
A  153 
A  154 
A  155 
A  156 
A  157 
A  158 
A  159 
A  140 
A  161 
A  162 
A  163 
A  164 
A  165 
A  166 
A  167 
A  168 
A  169 
A  170 
A  171 
A  172 
A  173 
A  174 
A  1*5 
A  176 
A  177 
A  178 
A  179 
A  180 
A  181 
A  182 
A  183 
A  184 
A  185 
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PRINT  81.  T.X.Y.Z.ALTfGAMMA.VMAS.PSI.V.QD.JtHACH.CD.ANR.PHltTHEtALP  A  186 

1>8ET,PS.ALPT,X£,YE,2E.FD,FN*GFZ.CLAT.DL0N.ABETA,GCIR,AIZ.AIX  A  18? 

LLLaLLL.l  A  188 

GO  TO  9  A  189 

8  D2X«FX/VMAS-2«*nA-C*WEY*B*KEZ*GX  A  190 

02Y»FY/VMAS-2.*nB-A*WEZ*C*WEY*GY  A  191 

02Z*FZ/VMAS-2.*nC-B**EX*A«¥EY*GZ  A  192 

02TH»TY/AIX««X2*0Z2*AIZ#BZ2*0X2/AIX-DBEY2  A  193 

02PH- {-TX/AlX*AIZ#WZ2*0Y2/AIX“BY2*0Z2*0PH*DTH*SlN{TH) ♦0*EX2) /COS (T  A  19* 

1H)  A  195 

02PS*TZ/AIZ-AlX«WY2*0X2/AIZ*AlX*WX2«0Y2/AlZ*D2PH*SlN(TH> *DPH*0TH*C  A  196 

10S<TH}»0tf£Z2  A  197 

CALL  COSF  (CD* AMACK, ANR)  A  198 

PER* 1 • ♦ • 0375»ALPT ♦ , 0  0 1 56*ALPT#*2 .  A  199 

C0*CD*P£R  A  200 

9  Xl«X^0X*H/2.  A  201 

DXl«0X*02X«H/2.  A  202 

Yl«Y*nY*H/2.  A  203 

Dy1»DY*02Y«H/2.  A  20* 

Zl*Z*0Z*H/2.  A  205 

0Zl-0Z*02Z*M/2.  A  206 

THl«TH*DTH«H/2.  A  207 

0TM1«0Tm*02TH«H/2.  A  208 

pmi«pm*0PH*H/2.  A  209 

0PMl«0PH»D2PH»M/2.  A  210 

PSl»PS*0PS»M/2.  A  211 

0P51-0PS*02PS*H/2.  A  212 

X5-X1  A  213 

0X5»0X1  A  21* 

Y5»Y1  A  215 

DY5»0Y1  A  216 

Z5*Z1  A  217 

0Z5-0Z1  A  218 

TH5«TH1  A  219 

0TH5-0TH1  A  220 

PH5«PH1  A  221 

0PH5-0PM1  A  222 

PS5-PS1  A  223 

0PS5«0PSl  A  22* 

01*1.800  A  225 

GO  TO  12  A  2?6 

10  D2X1«FX/VHAS-2.*0A-C*¥EY4B*MEZ*GX  A  227 

02Y1«FY/VMAS-2.*0B-A*MEZ*C«WcX*GY  A  228 

02ZI»FZ/v«AS-2,*0C-B»wEX*A*BEY*GZ  A  229 

02Thi«TY/AIX-WX2*CZ2*AIZ*MZ2*0X2/AIX-0KEY2  A  230 

02PhI«{«TX/A1X*aIZ*WZ2»0Y2/AIX-XY2*OZ2*0PM5*DTM5*STM5*OWEX2>/CTH8  A  231 

D2PSl«TZ/AIZ-AIx*WY2*0X2/AlZ*AIX*¥X2*0Y2/AIZ*D2PHl*STH5*DPH5*0TH5*  A  232 

1CTH5-0WEZ?  A  233 

X2»x*0Xl*H/2,  A  23* 

DX2»DX*02X1«h/2.  A  235 

Y2«Y*0Yl*M/2.  A  236 

0Y2»DY*02Yl*H/2.  A  237 

Z2*Z*nZl*H/2.  A  238 

DZ2*0Z«02Z1*h/2.  A  239 

TH2»TN*0THl*H/2.  A  2*0 

DTH2«0Tp*02Tm1«w/2,  A  2*1 

PH2«PM*nPHi*M/2.  A  2*2 

0PM2«0PH*02PHl»H/2.  A  2*3 

PS2*PS*0PS*  *N/2.  A  ?** 

0PS2*0PS*O2PS1*h/2.  A  2*5 

X5»X2  A  2*6 

0X5.0X2  A  2*7 
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11 


12 

13 


Y5«Y2 

0Y5*0Y2 

Z5«Z? 

DZ5*DZ2 

TH5*TH? 

0TH?«0Th3 

PH5*PH2 

0PH5*0PH? 

PS5»PS2 

0PS5»0PS? 

01*2.000 

GO  TO  13 

02X2*FX/VMAS-2 . •DA-C»WEY*8«*EZ*GX 
02Y2*FY/VKAS-2,.*DB-A*W£Z*C«*EX*GY 
02Z2*FZ/V«AS-2 . «0C-6*WE**A»lfEY*GZ 
0?T«2*TY/AlX-MX2*0Z2*AlZ«rt2»0X2/AIX-0KEY2 

02PM2*(.TX/AIX*4lZ*«Z2*OY2/ArX-WY2*OZ2*OPH5*DTHS*STHS*OMEX2)/CTHS 

02PS2*TZ/A1Z-AIx*MY?*OX?/.aI2*AIX«KX2*0Y2/AIZ*02PM2*STHS*DPH5*0TM5* 

1CTHS-0KEZ2 

X3*X»0X2*M/2. 

0X3*DX*02X2*H/2. 

Y3*Y*DY2*M/2. 

0Y3*0Y*02Y2*H/2. 

Z3*Z*nZ2*H/2. 

D73*DZ*02Z2*H/2. 

TH3*TH*0TH?*H/2. 

0TM3«0TM*0?Th2#h/2. 

pm3»ph*oph2*h/2 . 

0PH3«0Pm*0?Ph2*h/2 . 

PS3»PS*OPS2*m/2. 

DPS3«0PS*02PS2*H/2. 

X5«X3 

DX5*DX3 

Y5*Y3 

0Y5.*0Y3 

Z5*23 

0Z5*L'Z3 

TH5*1>’3 

0TH5*0r«3 

PMS*PH3 

0PH5*0P»3 

PS5*PS3 

0PS5*0PS3 

01*3.000 

T*T*DT 

RaSQRT < f RXH*X5) **2* (RYH*Y5) #*2* (RZH*Z5) **2; 

ZE* (RZH*ZSI ♦S6LAT-CGLAT* f fRXH*X5) •SAZ- fRYH* Y5) *CAZ> 

AL*ASIN(ZE/R) 

CALaCOS (AL1 
SAL*SIN(AU 
CL*T*AL*180./PI 

PE*REQ/ ( 1 .000*PE«SAL**2) ••♦500 

alt*r-re 

ye* (RXH.X6) *CAZ* <RYM*Y5)  *SAZ 

XE* (RZH*Z5) *CGLAT*SGLAT* < (RXH*X5>  «$AZ- IRYM* Y5) «CAZ) 

CALL  LONG  fOLG.XEtYE.PU 
CL6*C0SC0LG) 

SLG*SINfOLG) 

0L0N*0LGG*0LG*1A0./PI  _  .  . 

GCIR*( (RE*REl )/2.)*AC0S(SAL*5lN(CLATl) ♦CAL#C0SfCLATl )»C0S(AL*CLAT1 

1)1 

G11*-SGLAT*CAL*CLG*CGLAT*SAL 


A  248 
A  249 

A  250 
A  251 
A  252 
A  253 
A  254 
A  255 
A  256 
A  257 
A  258 
A  259 
A  260 
A  261 
A  262 
A  263 
A  264 
A  265 
A  266 
A  267 
A  268 
A  269 
A  270 
A  271 
A  272 
A  273 
A  274 
A  275 
A  276 
A  277 
A  278 
A  279 
A  260 
A  281 
A  282 
A  263 
A  284 
A  285 
A  286 
A  287 
A  266 
A  289 
A  290 
A  291 
A  292 
A  293 
A  294 
A  295 
A  296 
A  297 
A  298 
A  299 
A  300 
A  301 
A  302 
A  303 
A  304 
A  305 
A  306 
A  307 
A  308 
A  309 
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612«SGLAT»S*L*Cl6«C6LAT»CAL 

G21»-CA|_*SL6 

G22«S*L*SLS 

63 1 ■-CSL*T*C*L«CLG-S6L*T«S*L 
632»*C6L*T*S*L»fiL6-S6L*T»C*L 
P12«1.-3.*S*L*»? 

P13*J.*S*L-5.»SAL**3 

P1*.3.-30.*SAL*«2O5.*SAlp** 

PlS.SINl2.»AU 
PS *CAL* 11.-5 . »S*L*«2 ) 

P?«SAL*CAL* l -3 . *?.«SAL**2  > 

G*».U/R*«2* C 1 .*8 J* (R£8/R) ••2«P12*4./5.*SH« (REQ/R) ••S*.'13*3K/6.« {»£ 
ia/Ri****Pi4j 

Si-*U/R*#2*  J*{R£Q/Rl**?*®15«3»/5.*RM*  tR£Q/RJ4«34P6*2./3.*8X*fR£0 

6x*SR* (611*5A2»621*C*Z) •  6LM612*S*Z*S22*CAZ) 
GY»R»*(S21*SAZ-G11*CAZ>*6LM6Z2*SAZ-S12*CAZJ 
SZ«fiR»S31*6L*S3? 

CPS5^C0S{PS51 

SPS«5«SlNtPS5) 

STHS.SINIThS) 

CTH5*C0SITm5} 

SPH6«$I*{Pm5i 
CPHS*C0S(PK5l 
A 1 1 •CP55*CTHS 

A1 2«S®S5»CPH5-STH5*SPM5«C*>55 
*1 3«-CPS5»STh5*CPH5-SP5S*RPH5 
A21.-SP55*CTk5 
A22.S?S5*STt«5«SPK5*CP55PC*H5 

A23»SPS5*STH5*CPH5-SPH5*CPS5 

A31«RTH5 

A32*CTH5»SP«5 

A33«CTH5»CPH5 

a*hey* {Q2H-Z51 - (RYH-Y5) **EZ 
B.tfEZ*  t RXK-X5) - f RZM*Z5 J  *¥E* 

C«tfEX«  tRTH*T51 - c RXM*X51 «*¥£? 

OA*OZ5*«EY-Ot5*¥EZ 

OR»nx5»lfEZ-OZ5*«EX 

oc*nY5«t»£x-Dx5«i^;r 

8ll«CTH? 

*'12»-STM5*SPhS 
u  j  3«-ST  H5«C?*5 
*21»0. 

822.CPH5 

8?3«— SPmS 

331.STH5 

932»CTM5»SPH5 

833»rCTM*;«CP«5 

D8U»-DT*<?*ST*S 

03 1  2*-Dth5»CT«5»SPiS-T>PH5»CPk5*STh3 
0R13*-0Tk5*CTH5»CPms»0PH8«^m5*ST»«5 
0*21.0. 

D322.-3Pm««SPm5  ^7 

0a23.-0P«5*CPK5 

0331«DTm5*CTk5 

D332«0Pw5*CPh5*CT*5-DTh5«STK5*SP**5 
DA33«-0THS»STh5*CPm5-DPX5«SPH5*CT»<5 
WEX2.Sll*ifEX.9i2*irEY*813*KEZ 
*EY?.32 1  •wEX*e2p*t(EY*SP3»wEZ 
*EZ?^3l»ifEX-83P*»fEV*a334»EZ 
DWEX2.0R  1 1  •*EX-S912«*EY*0R1 3*lr£Z 
0*Er?.Ds2l  •«EX*is322*Arr*J>R23*«EZ 


A  330 
A  311 
A  312 
A  313 
A  314 
A  335 
A  316 
A  31? 
A  31 A 
A  319 
A  326 
A  321 
A  322 
A  323 
A  324 
A  325 
A  326 
A  327 
A  3?8 
A  329 
A  330 
A  331 
A  332 
A  333 
A  334 
A  335 
A  336 
A  337 
A  338 
A  339 
A  340 
A  341 
A  342 
A  3*3 
A  3*4 
A  345 
A  3*6 
A  347 
A  3*8 
A  3*9 
A  350 
A  351 
A  352 
A  353 
A  35* 
A  355 
A  356 
A  357 
A  358 
A  359 

*  3*0 
A  361 
A  36? 

*  363 
A  36* 
A  365 
A  366 
A  367 
A  368 
A  369 
A  37g 
A  371 
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n  o 


DWE7?*DR31*WEX*DS3?«WEY*0B33»#EZ  A  372 
WX2«WFX2-DPN5*CTH5  A  373 
WY2«WEY2*DTH5  A  37* 
WZ2*WFZ?*nPSS-D0H5*STH5  A  375 
0X2-WX2  A  376 
0Y2*wY2  A  377 
0Z2*WEZ?-DPH5*STH5  A  378 
VPX-0X5  A  379 
v/RY«DY5  A  380 
VRZ«nZ5  A  381 
VR*SQRT  (VRX**?*VRY#-'24VRZ**2)  A  382 
VsSORT(0X5**2*DyS<*#2*D7S##2)  A  383 
CALL  RHOr  (HHO» ALT )  A  38* 
CALL  VISCO  (VIS, ALT)  A  385 
gn*,«;*RHO*V**2  A  386 
CALL  VSD  (VS, ALT)  A  387 
AMACH«V/VS  A  388 
ANR»RHO*V«SLANT/VIS  A  389 
VRX*»A11*VRX*A12*VRY*A13*VRZ  A  390 
VRYAw£2 1 *VRX* A22*VRY*A23*VRZ  A  391 
VRZ*>A31*VRX«A32#VRY«A33*VRZ  A  392 
WX3«-0Ph5*CTH5*CPS5*0TH5*SPS5  A  393 
WY3«0PH5*CTH5*SpSS*DTH5»CPS5  A  39* 
WZ3»OPS5-nPH5*STH5  A  395 
CALL  ATTK  ( VRZ*. VRYA.ALP.PI)  A  396 
CALL  SSLP  ( VRZ4, VRX*,BET ,Pl )  A  397 
ALPTbALS(ALP)  A  398 
ALPT»ALPT*180./pI  A  399 
ALP«ALP*180./PI  A  *00 
8ET«BET*180./PI  A  *01 
ALP1»AB«!F(ALP)  A  *02 
BETl *A8SF (BET )  A  *03 
CALL  ACMO  (CMQ.ftMACH)  A  *0* 
CALL  ACNA  (CNA.AMACH)  A  *05 
CNA  PER  DEGREE  A  *06 
CMQ  PER  RAOIAN  A  *07 
CNRaCMO  A  *08 
CMAa-SM*CNA/DIA  A  *09 
CM«CMA*ALP*CM0*WX3*DIA/(2,*V)  A  *10 
TM»CM*QD*AREA«DI A  A  *11 
FN*CNA*QD*AR£A*ALP  A  *12 
CNBa-CMA  A  *13 
CY0«-CNA  A  *1* 
CN»CNB*BET*CNR«WY3*0IA/(2,*V)  A  *15 
TN»CN*Qn*AREA*DIA  A  *16 
FY1»CYB*Q0*AREA*8ET  A  *17 
TL«0.  A  *18 
F0«CD*QD*AREA  A  *19 
FX3«FY1  A  *20 
Fy3»-FN  A  *?1 
FZ3»-FD  A  *22 
GFX»FX3/(VMAS*GRAV)  A  *23 
GFY»(FY3-VMAS«*GrAV*SJN(PH)  )/(VMAS*GRAV)  A  *2* 
GFZ«(FZ3*VMAS*GRAV»C0S(PH) ) MVMAE^GRAV)  A  *25 
ABETA«VMAS<*3RAV/(CD*AREA)  A  *26 
FX*AU*FX3*A21«*FY3*A31*F23  A  *27 
FYsA12«FX3^A?2*FY3*A32*FZ3  A  *28 
FZ»*13*FX3*A-  3«FY3*A33*F73  A  *29 
TX3«TM  A  *30 
TY3»Tm  A  *31 
TZ3-TL  A  *32 
TX«TX3»CPS5-TY3*SPS5  A  *33 
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F 

-  -  -  — 

-  ■  “  WU.»XV>lWj?  L  .  -■-■**111  IIP  III  -■  I 

-  - 

: 

i 

; 

f  i 

JT 

it 

? 

• 

\ 

5 

t 

j 

TY«TX3*SPS5*TY3*CPS5 

A 

434 

T7-TZ3 

A 

435 

IF  ’4.000-01)  15*18.14 

A\ 

436 

u 

IF  (1. 000-01)  14*10*15 

A 

437 

*  15 

GO  TO  69 

A 

438 

16 

IF  (2.000-01)  17.11*15 

A 

439 

IT 

D2X3*FX/VMAS-2.*DA-C*WEY*a*¥EZ*GX 

A 

440 

* 

D2Y3eFY/VMAS-2.*D3-A*WEZ*C*«EX*GY 

A 

441 

02Z3aFZ/VMAS-2.*DC-8*WEX*A*«EY*GZ 

A 

442 

02TH3aTY/AIX-WX2*0Z2*AIZ*«2*0X2/AlX-DWEY2 

A 

443 

D2PH3*(-TX/A!X*AIZ*WZ2»0Y2/AIX-WY2*0Z2*DPH5*DTH5*STH5*DWEX2)/CTH5 

A 

444 

D2PS3*TZ/-IZ-AIX*WY2*0X2/AIZ*AIX*WX2«0Y2/AIZ*D2PH3*STH5*0PH5*DTH5* 

A 

445 

1CTH5-0WF.Z2 

A 

446 

OELXaH/6.* (DX*2 ,*DX1 *2»*DX2*DX2) 

A 

447 

0ELDXaH/6.*(02X*2. *02X1 *2. *02X2*02X3) 

A 

448 

; 

DELYaH/6.* (DY*2.*0Y1 *2.*0Y2*DY3) 

A 

449 

0EL0y*H/6.* (D2Y*2,*02Y1*2«  *02Y2*02Y3) 

A 

450 

DELZ«H/6.*(0Z*2.*0Z1*2.*0Z2*0?3) 

A 

451 

; 

DELOZ-H/6,*(D2Z«2.*D2Z1*2,*02Z2*02Z3) 

A 

452 

DELTH«H/6.*(0TH*2.*0TH1*2,*0TM2*0TH3) 

A 

453 

DELOTH*H/6,*(02TH*2.*02TH1*2.*D2TH2*02TH3) 

A 

454 

* 

- 

0ELPH«K/6.*(0PH*2.*0PH1*2.*0PH2*DPH3) 

A 

455 

0EU0PH-h/6.*(D2PH*2.*02PH1*2.*02PH2*02PH3) 

A 

456 

DELP5-H/6.*(DPS*2.*0PS1*2.*0P$2*DPS3> 

A 

457 

0EI-0PS»H/6.*(D2pS.2.*02PSI*2.*D2PS2*D2PS3) 

A 

458 

• 

X»X*OELX 

A 

459 

0).*0X*0EL0X 

A 

460 

Y»'Y*OELY 

A 

461 

’ 

dy«oy*oeloy 

A 

462 

z*z*oelz 

A 

463 

X 

dz*dz*deloz 

A 

464 

TH«TH*OELTH 

A 

465 

■ 

DTH»OTH*OELOTH 

A 

466 

ph*ph*oelph 

A 

467 

! 

oph«dph*oeloph 

A 

468 

PS-PS*DELPS 

A 

469 

f  | 

0PS«0PS*0EL0PS 

A 

470 

1  1 

X5*X 

A 

471 

\  \ 

DXSaDX 

A 

472 

i  * 

Y5»Y 

A 

473 

r 

D~>  5*DY 

A 

474 

I 

Z5«Z 

A 

475 

i  i 

DZSaOZ 

A 

476 

f  i 

TH5aTH 

A 

477 

I 

DTH5«0TH 

A 

478 

PH5*PH 

A 

479 

i 

OPHSaOPH 

A 

480 

E 

PS5aPS 

A 

4R1 

§E 

OPSBaOPS 

A 

482 

§? 

91*4.000 

A 

483 

I 

GO  TO  13 

A 

484 

-* 

16 

THEaTH*180./PI 

A 

485 

PHlaPH*18C,/PI 

A 

486 

£ 

GAHMAaATAN (0Z5/DY5) *180./PI 

A 

487 

B 

PS JaATAN (QX5/DY5) *180 ./PI 

A 

488 

S  £ 

8TLaBT-H/4, 

A 

489 

'  i 

BTHaBT*H/6. 

A 

490 

% 

IF  (ALT. LE. 0.0)  GO  TO  20 

A 

491 

S. 

IF  (T.ST.RTL.ANn.T.LT.BTH)  GO  TO  19 

A 

492 

1 

GO  TO  8 

A 

493 

i 

19 

BT»8T* .05 

A 

494 

C 

A 

495 

8 

99 

It 

* 

i 

A  496 


WL-0.0  A  497 

00  49  LENa  1  »N0FX  A  498 

WL«WL*SCL(LEN)  A  499 

IF  (LEN.LT.NOFX)  GO  TO  21  A  500 

WL-.25  A  501 

NALaNALC  A  502 

00  22  Nal.NAL  A  503 

ANaN  A  504 

IF  (AN*DEL<1.1)-S00T(LEN5)  22t22*23  A  505 

CONTINUE  A  506 

PRINT  84  A  507 

GO  TO  69  A  508 

KaN  A  509 

IF  (LEN.GE.N0FX)  TMELT«64d0.  A  510 

TR«TENPE*.58*(Tw(K,LEN)-TEMPE)*.19MTAtf-TEMPE>  *  511 

CALL  TEMPA  (TEMP. ALT)  A  512 

CALL  TCAT  (TCA.TR)  A  513 

CALL  VISCOT  (VISE.TEMPE)  A  514 

CALL  VISCOT  (VISG.TR)  A  515 

CALL  PRESS  (PRE.ALT)  A  516 

AMaAMACH  A  517 

PBLaPRE* ( 1 .♦2.8«AM«*2* (SIN(THETC) ) #*2* (2.5^8.#AM»SIN (~HETC» ) / ( 1 .*  1  A  518 

16.*AM«SIN(THETC) J )  A  519 

P»P3L/21 16.2  A  520 

IF  (JIND(LEN).EQ.OJ  GO  TO  24  A  521 

•••••••••••*  a  522 

CALL  GCP8  (uCP.P)  A  523 

IF  (LEN.GE.NOFX)  3CPa,4  A  524 

a  525 

CALL  GASTC  (GTC.GCP)  A  526 

PRaGCP*GVI5/GTC  A  527 

GO  TO  25  4  528 

PR«CPA*VISG/TCA  A  529 

HNRFaRHOE*VE»WL/V!SE*G  A  530 

RH0S»PBE/(l7l6.«TR)  A  531 

ANRRaRHOS*VE*WL/VISG*G  A  532 

IF  (HNRP.GT.RN(LEN) )  GO  TO  26  A  533 

REaPRa«.5  A  534 

HFa.575«TCA/WL«PR**.33«HNRE*».5*<VlSG/VISE*RH0S/RH0E)*».5  A  535 

GO  TO  27  A  536 

RFaPR»«.33  A  537 

HE».0296*TC4/XL*PR*#.33#MNRE#*.8# (VISG/VISE) ##.2# (RH05/RM0E) **.8  A  538 

0TRavE**2/ (2,*G*AJ*CPA)  A  539 

TAtoaTEMP£.RF»OTR  A  540 

OOOTC(LEN)aHE*(TAW-TW(K»LEN) )  A  541 

IF  (INO(K.LEN) .EQ.1i  GO  TO  33  A  542 

IF  (Th(K.LEN)-TmELT)  29.33.33  A  543 

TWKaT*VK.LEN)  A  544 

IF  (LEN.GE.NOFX)  GO  TO  30  A  345 

CALL  VTCPC  (VT.TVK)  A  546 

CALL  VCPPC  (VC.TWK)  A  547 

VRHO (K.LEN) aRHOV  A  548 

GO  TO  31  A  549 

CALL  TCATJ  (VT.TWX)  A  550 

CALL  CPATJ  (VC.TWK)  A  551 

VRHO (K.LEN) a AT JRHO  A  552 

VTC (K »LEN) aVT  A  553 

VCP(K.LFN) a VC  A  554 

IF  (K.EQ.l)  GO  TO  36  A  555 

Pl«VTC(K-l.LEN)«(TW(K-l.LEN)-TW(K.LEN> )/(VRHO(K  l.LEN)*VCP(X-l»LFN  A  556 

1 ) *0E ( 5-1 *LEN) ••2)  A  557 


32 

33 


3* 

35 

36 


37 

38 


39 

*0 

*1 

*2 


43 


44 

45 

46 


P2«VTC(K«LEN)*(TW(K,LEN)-TW(K*1.LEN) )/ (VRHO (K.LEN) *VCP (K.LEN) »DE (K 
l.LEN)*#2) 

TW<K,LEN)«TK(K,LEN)*DELT»(P1-P2) 

IF  (TW(K.LEN) -TMELT)  48*32*32 
TW(K,LEN)«THELT-1. 

60  TO  48 
TWK«TW (K.LEN) 

IF  (LEN.GE.NOFX)  GO  TO  34 
CALL  CTCPC  (CT.TWK) 

CALL  CCPP  (CC.TWK) 

VRHO (K.LEN) *RHOC 
GO  TO  35 

CALL  TCATJ  (CT.TWK) 

CALL  CPATJ  (CC.TWK) 

VRHO (K.LEN) «ATjRHO 
VCP (K.LEN) «CC 
VTC (K.LEN) »CT 
INO (K.LEN) al 

33-DELT/ ( VCP (K.LEN) *VRHO (K.LEN) #0E (K.LEN) ) «*2. 

P4«HE« (TAW-Tw (K.LEN) ) 

P5*vTC (K.LEN) *  (TW (K.LEN) -TW  <K* 1 .LEN) > /OE (K.LEN) 
Hp«CPA*riNPE*RF»VE**2/ (2.»G«AJ) 

QCONO (LEN) «-VTC (K.LEN) • t TW (K.LEN) -TW (K*l ,L£N) ) /DE (K.LEN) 

TW (K.LEN) -TW (K.LEN) *P3* (P4-P5) 

IF  (TW (K.LEN) -TMELT)  37.38.38 

JINO(LEN)«o 

WT(L£N)«0, 

GO  TO  48 
TW (K.LEN) -TMELT 

IN0(K,LEN)«1 
JINO (LEN) »1 
CALL  GCP8  (GCP.P) 

IF  (LEN.GE.NOFX)  QCP*.4 
QDOTR-0. 

OTOT (LEN) «QDOTC (LEN) *QCOND (LEN) 

IF  (QTOT(LEN)  .LE.O.O)  GO  T,;  48 
IF  (LEN.LT.NOFX)  GO  TO  39 

delh*hr-cpa*tmelt 

QST4R»2000.*2.»OELH 
COOT.QTOT  (LEN)  /  (QSTAR*AT-;rxHO) 

GO  TO  4' 

IF  (HNRE-RN(LEN) )  40.40.41 

0MD« (QTOT (LEN) ) / ( AKL1 *AKL2* (MR-CPA«TMELT) > 

GO  TO  42 

OMO» (QTOT (LEN) ) / ( AKTl*AKT2« (HR-CPA*TMELT) ) 

PXE«11.05E*04/TMELT 

TMDOT«O«D*(l.*2.64E*09/(P8L*».674E**PXE) ) 

COOT-TMOOT/RHOV 

SHOT (LEN) «SOOT (LEN) *CDOT»DELT 

JINO(LEN)«l 

00  44  N-l.NAL 

BN«N 

IF  (8N*0EL (1.1) -SOOT (LEN) )  44.44.45 
CONTINUE 
PRINT  84 
GO  TO  69 

HE (N, LEN) »0EL (NtLEN) - (SOOT (LEN) - (BN-1 . ) *OEL ( 1 .LEN) ) 

PCM»VCP (N.LEN) *VRHO (N. LEN) *0E (N.LEN) «#2/ (OELT#VTC (N.lEN) 5 
IF  (PCM-2.)  46.47.47 
OE (N.LEN) »n  * 

CN«N 

8DOT(LEN)»CN*TOL 


A  558 
A  559 
A  560 
A  561 
A  562 
A  563 
A  564 
A  565 
A  566 
A  567 
A  568 
A  569 
A  570 
A  571 
A  572 
A  573 
A  574 
A  575 
A  576 
A  577 
A  578 
A579 
A  5«0 
A  581 
A  582 
A  583 
A  584 
A  585 
A  586 
A  587 
A  588 
A  589 
A  590 
A  591 
A  592 
A  593 
A  594 
A  595 
A  596 
A  597 
A  598 
A  599 
A  600 
A  601 
A  602 
A  603 
A  604 
A  605 
A  606 
A  607 
A  608 
A  609 
A  610 
A  611 
A  612 
A  613 
A  614 
A  615 
A  616 
A  617 
A  61  6 
A  619 


101 


» 


IND(N*1,LEN)«1 

IND(N,LEN).l 

kT (LEN) «PI*  (RL (lEN) *RS (LEN) -2.*SD0T (LEN) ) «COOT*RHOC*SOL (LEN) 

IF  {N.E0.1J  60  TO  48 
0E(N-1»LEN)»0. 

K*K*  1 

IF  (K.LE.NAL)  GO  TO  28 

TW<K,LEN)«TK(K,LEN)*DELT*(VTC(K-1.LEN)/(VRH0{K-1,LEN)*VCP(K-1,LEN) 

1 *0E (X-l .LEN) **2) * (Tk  (K-l «L£N) -Tk (K.LEN) ) -VTC (K.LEN) / (VRHO(K.LEN) *V 
2CP  <K. LEN) *OEL (K.LEN) ••Z) • (Tk (K.LEN) -Tk (K* 1 .LEN) ) ) 

IF  (TW(K.LEN) .Gt.1060.)  Tk <K,L£N) «1060. 

K*K*1 

TW(K.LEN)«Tk{K.L£N)*DELT*(VTC{K-l.LEN)/(VRHO(K-l,LEN)*VCP(K-l.LEN) 
1*0EL(K-1.LEN)**2)#(TW(K-1,LEN)-TW (K.LEN))) 

IF  (Tk(K.LEN) .GT.860.)  Tk (K.LEN) «86C. 

NALaNALI 

CONTINUE 

LENbNOFXP 
00  50  N.l.NALS 
AN«N 

IF  (AN»0EL(1.N0FXP)-S00T(LEN) )  50.50.51 

continue 

PRINT  85 
GO  TO  6G 
KS»N 

TMELT«6*00. 

STR»TEMP*.S8«(Tw<KS.LEN)-TEMp)*.19*(STAk-TEMP) 

CALL  TCAT  (STCA.STR) 

CALL  VI5C0T  (SVISG.STR) 

IF  (JINO(LEN) .EO.O)  GO  TO  52 

CALL  GCPB  (GCP.P) 

IF  (LEN.GE.NOFX)  GCP-.4 

••••••••••••••••••••••«•••••••••  .>***••••*•••*«***•***«•**•••**•• 

CALL  GASTC  (GTC.GCP) 

SPR»GCP*GV I S/GT C 
GO  TO  53 

SPR-CPA«SVISG/STCA 

SRF«SPR«*.5 

CALL  PRESS (PRE. ALT) 

PTS»PRE# ( O - "*AmACH«*2)#*3.5)*( (6,0/(7.0»AkACH**2-1.0) >##2.5) 
RHOkS-PTS  i j.*STR) 

DV»1./SN'  ■'  1 'O.*(PTS-PRE)/RH0kS) 

GTAB« ( ( t .  ,+S.C/ ( AMACH**2) ) * ( 1 . 0-1 . 0/ ( 1 . A*AMACH**2) ) ) **0.25 
SGDOTC-0.<^5--i3*GTAB*S0RT(RH0kS«*SVISG*DV»G  )  •  (V*»2)  / (2.*G*AJ) 
SOTR«V'**2/  (2.#G»AJ*CPA) 

STAk«TEMP*SRF*SnTR 

stg«temp*sotr 

STREC«1EMP*SRF«(STG-TEMP) 

SHE»SQOOTC/(STG-Tk(KS.LEN) ) 

IF  (IND(KS,LEN).£0.1)  GO  TO  57 
IF  (Tk(KS.LEN) .GE.TNELT)  GO  TO  57 
STkKaTk (KS.LFN) 

CALL  TCATJ  (SVT.STkK) 

CALL  CPATJ  (SVC.STkK) 

SARH0(K5)»ATJRH0 
SATC(KS)»SVT 
8ACP (KS ) »SVC 
IF  (KS.EO.l)  GO  TO  58 

Tk(KS«LEN)»TW{KS.LEN)*DELT*(SATC(KS-l>/(SARHO(KS-l)#SACP(KS-l)«DE( 
lKS-l.LEN)*«2)*(Tk(KS-l.LEN)-Tk(KS.LEN) )-S*TC(KS)/(SARHO(KS)*SACP(K 


A  620 
A  621 
A  622 
A  623 
A  624 
A  625 
A  626 
A  627 
A  628 
A  629 
A  630 
A  631 
A  632 
A  633 
A  634 
A  635 
A  636 
A  637 
A  638 
A  639 
A  640 
A  641 
A  642 
A  643 
A  644 
A  645 
A  646 
A  647 
A  648 
A  649 
A  650 
A  651 
A  652 
A  653 
A  654 
A  655 
A  656 
A  657 
A  658 
A  659 
A  660 
A  661 
A  662 
A  663 
A  664 
A  665 
A  666 
A  667 
A  668 
A  669 
A  670 
A  671 
A  672 
A  673 
A  674 
A  675 
A  676 
A  677 
A  678 
A  679 
A  680 
A  681 


¥ 


i 


57 


5R 


59 


60 

61 

62 


63 


*.4 

65 


C 


2S) *OE (KS.LEN)  ••?) • (TW (KS.LEN) “TW (KS*1 »LEN) ) ) 

IF  (TW(KS.LEN) .LE.TMELT)  60  TO  65 
TW (KS.LEN) aTKELT-1. 

60  TC  65 
STHKaTN (KS.LEN) 

CALL  TCATJ  (SCT.STKK) 

CALL  CPATJ  (SCC.STWK) 

SARHO (KS)*ATjRHO 
SATC(KS)«SCT 
SACP (KS) «SCC 
IND{KS»LE*0»1 

TW (KS.LEN) aTfc (KS.LEN) *OELT/ (SACP (KS) *SARHO (KSi *DE (KS.LEN) ) * (SHE* (S 
1T0  -TW {KS.LEN) ) -SATC  (KS) /DE  (KS. LEN)  MT8(KS. LEN) -T«(KS*1. LEN) )  )*2. 
IF  (TW(KS.LEN) .GE.THELT)  60  TO  59 
JINO(LEN>«0 
GO  TO  65 
TW(*S*LEN)«TMELT 
HR«CPA*TE^-v*RF*V**2/  (2»*6*AJ) 

SOCONO—SATC  (KS) • (T* (KS.LEN)  -T*  (KS*1  .LEM) )  /DE  (KS.LEN) 

SQTOT-SOOOTC.SQCONO 

IF  (SOTOT.LE.o.fl)  GO  TO  65 

IND (KS.LEN) «1 

JIN0(LEN)«1 

oelh«hr-cpa«tmelt 

OSTARS*2000.*2.*DELH 
COOTSaSOTOT/ (OSTARS*ATJRHO> 

SOOT (LEN) -SOOT (LEN) .CDOTS^OELT 
SS-CD0TS 

RX»RX.1./(1.-SIN(TMETC)>*(SIM(THETC)«SS-CD0T) 

SX».01745*RX*PMI 
IF  (SX.LE.O.O)  SX«. 01943 
IF  (RX.LE.O.O)  RX«. 01387 
SNRaRX 

00  60  Nal.NALS 
BNaN 

IF  (8N*0EL(1.N0FXP)-S00T(LEN))  60.60.61 
CONTINUE 
PRINT  84 
GO  TO  69 

OC (N.LEN) »OEL (N.LEN) - (SOOT (LEN) - (BN-l . >  *OEL ( 1 .NOFXP) ) 

PCH»SACP (N) «SARHO (N) *DE (N.LEN) **2/ (DELT*S*TC (N) ) 

IF  (PCM-2.)  62.63.63 
DE(N,LEN)»0. 

SOOT (LEN) .BN^STOL 
IND (N.LEN) al 
IND (N*l .LEN) al 
IF  (N.EG.l)  00  TO  65 
NMaN-1 

00  44  Ial.NM 
OE(I.LEN)aO. 

KSaKS*l 

IF  (KS.LE.NALS)  GO  TO  56 

Tw (KS.LEN) aTW (KS.LEN) *DELT* (SATC CKS-1 ) / (SARHO (KS-1 ) *SACP (KS-1 ) *0F ( 
1KS-1 .LEN) **2) • (TK (KS-1 .LEN) -T M (KS.LEN) ) -SATC (KS) / (SARHO (KS) *SACP  (K 
pS)*0E(KS.LEN)**?)*(TW(KS.LEN)-T»#(KS.1.LEN)  ) ) 

IF  (TW (KS.LEN) .ST. 1060. )  T* (KS.LEN) al 060 . 

KSaKS.l 

TW (KS.LEN) aTK (KS.LEN) .CELT* (SATC (KS-) 1 / (SARHO (KS-1 ) *SACP (KS-1 >  *DE ( 
1KS-1. LEN) **2)*(TW (KS-1 .LFN)-TK (KS.LEN) ) ) 

IF  (TW(KS.LEN) .GT.860.)  Tk (KS.LEN) «860 . 


A  682 
A  683 
A  684 
A  685 
A  686 
A  687 
A  688 
A  689 
A  690 
A  691 
A  692 
A  693 
A  694 
A  695 
A  696 
A  697 
A  698 
A  699 
A  700 
A  701 
A  702 
A  703 
A  704 
A  705 
A  706 
A  707 
A  708 
A  709 
A  710 
A  711 
A  712 
A  713 
A  714 
A  715 
A  716 
A  717 
A  718 
A  719 
A  720 
A  721 
A  722 
A  723 
A  724 
A  725 
*  726 
A  727 
A  728 
A  729 
A  730 
A  731 
A  732 
A  733 
A  734 
A  735 
A  736 
A  737 
A  738 
A  739 
A  740 
A  741 
A  742 
A  743 


TMELT*6760 
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c  *******«****•*••••*******•*#•*»••••*•••*•*♦#••*•****••*•♦*••♦•••**  A  74* 

TOtWT*0  >  A  7*5 

WOOTbO.O  A  7*6 

00  66  LL*1.N0FX  A  7*7 

T0TWT»T0TWT*PI*  (PS  (i.1.)  ♦HL (LL>  “SOOT  C |_|_ i )  *SOOT  (LL) *RHOV*SOL  (LL)  A  7*8 

66  wOOT.WDOT.wT (LL)  A  7*9 

SHTL*.0698#pHl*SNRI**2*SOOT(NOFXP}*ATJRHO  A  750 

WT(NOFXP)»2.*PI*SNR**2*COOTS<*ATjRHO  A  751 

TOTwT»TOTwT*SwT|  A  752 

WO0T«WD0T.WTtN0FXP)  A  753 

VMAS«8I-T0TWT/G  A  754 

RAD«RADI-SOOT(NOFXH)  A  755 

ALGTH»ALGTHl-SOOT(NOFXP)  A  756 

AREA»PI*RAD#*2  A  757 

SLANT«ALGTH/ (C08(THETC) -SNR/RAO* (1. -SIN (THETC)> )  A  758 

ATL«AT-0ELT/6.  A  759 

ATH«AT*0ELT/6.  A  760 

IF  (ALT. LE. 0.0)  GO  TO  68  A  761 

IF  (T.GT.ATL.AND.T.LT.ATH)  GO  TO  67  *  762 

GO  TO  8  A  763 

67  IF  (LLL.LT.8)  GO  TO  68  A  76* 

LLL*0  A  765 

PRINT  80  A  766 

68  PRINT  81*  T.X,Y,Z*ALT.GAMMA«VMAS.PSI.V*0D.AMACH.C0.ANR.PH1,THE»ALP  A  767 

1 .BET »PSi ALPT.XE.YE.ZE.FD.FN.6FZ.CLAT.DL0N. A8ETA.6CIR.AIZ.AIX  A  768 

PRINT  82.  CP.CAN.CPB.COF.ANRE  A  769 

PRINT  83.  RHO.VIS.AMACHE.TEMPE.RHOE.VE.VISG.PBE  A  770 

LLULLL*)  A  771 

AT»*T* 1 . 0  A  772 

PRINT  87.  (SOOT(I) ,I*1,N0FXP)  A  773 

PRINT  70.  mw<T.J).J»l.N0FX).I«l.NP2)  A  77* 

PRINT  86*  TW(NALCl.NOFX) ,TW<NALC2.NOFX>  A  775 

PRINT  89.  WDOT.TOTNT.RX  A  776 

PRINT  90.  QSTAR  A  777 

PRINT  9),  (TW(I.NOFXP) .I-1.NP2S)  A  778 

IF  (ALT)  69,69.8  A  779 

69  CONTINUE  *  780 

C  A  781 

T  A  782 

70  FORMAT  (30X.8F8.0)  A  783 

71  FORMAT  (1H1)  A  78* 

72  FOPMAT  (*JIN0*10I*)  A  785 

73  FORMAT  <*IND»8I4/(10X»8I*) )  A  786 

7*  FORMAT  (•VCP*8El2.3/(llX,8E12.3) )  A  787 

75  FORMAT  (#VTC*R£l 2,3/ (11X.8E12.3))  A  788 

76  FORMAT  (*VRH0*8E12.3/(11X.8E12.3) )  A  789 

77  FORMAT  (*0£L*8E 1 2,3/ ( 1 1X.8E12.3) )  A  790 

78  FORMAT  (*DE*«E12.3/(11X.8E12.3) )  A  791 

79  FORMAT  (•Tw*8El?.3/ (11X.8E12.3) )  A  792 

80  FORMAT  (12X.*HTIME,1*X.1HX.18X.1HY.15X.1HZ.16X,3HALT,13X.5HGAMMA,1  A  793 

11X.*HVMAS/30X,*hPSI.1SX,3HVEL»13X.2MQD»15X,5HAMACH,1ix.2HC0,1*X,3H  A  79* 

2REN/30X,3HPHI,16X.5HTHETA,UX«5HALPHA.12x»3HBET,13X»2HPS,1*X.*HALP  A  795 

3T/30X.2HXE.17X.?HYE.l*X.2MZE»15X»2HFD»l*X,2HFN,l*X,3hGFZ/30X,3HLAT  A  796 

*,16X,*HL0NG,I2X,*HBETA.13X.*HGCIR.12X.2HIZ.1*X,2HIX)  A  797 

81  FORMAT  (1H0»*X.7E17.7/(22X,6E17,7) )  A  798 

82  FORMAT  (10X.5E1O.3)  A  799 

83  FORMAT  (5X,8E’0.3>  A  800 

8*  FORMAT  (10M  OVERHEAT/)  A  801 

85  FORMAT  (•0VERHEATSTAG1*/)  A  802 

86  FORMAT  (*0VERHFATSTAG2*/)  *  803 

87  FORMAT  ( 10X#S0OT*9E12.3)  A  80* 

88  FORMAT  (86X.F8.ft)  A  805 
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89  FORMAT  (10x«*TLOSSRATE!$«Ell.**LBS/SECANDTOTA!.*rLOSSlS*tl3.**L8SAN  A  806 

10N0SERA0IS*Ell.**FT*l  A  807 

90  FORMAT  ( 10X#0STaRIS*E1 1 •♦•BTU/LB*)  A  808 

91  FORMAT  (30X.13FR.01  A  809 

92  FORMAT!*  SOOT*  10F8.S)  A81C 

END  A811 
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SUBROUTINE  ATTK  (VRZ4.VRr4.ALP, PI ) 
C  COMPUTES  ANGLE  OF  A'TACK.ALP* 

IF  (VRZ4)  9.5*1 

1  IF  (VRY4)  4.3,2 

2  ALPsATAN ( VRY4/VRZ4) 

RETURN 

3  ALP«0. 
return 

ALP«ATAN(VSY4/VRZ4) 

RETURN 

IF  (VRY4)  6.7,8 
ALP—PI/2. 

RETURN 
ALP«0. 

RETURN 

8  ALP*PI/>. 
return 

9  IF  (VRY4)  10,11,12 

10  ALP»-PI*ATAN(VRY4/VRZ4) 
return 

11  ALP«PI 
RETURN 

12  ALP»P I • AT AN ( VRY&/VRZ4 ) 

RETURN 

END 


8  1 
B  2 
B  3 
B  4 
8  5 

8  6 
8  7 

8  8 
B  9 
B  10 
B  11 
B  12 
B  13 
8  14 

B  15 
8  16 
8  17 

B  18 
B  19 
B  20 
B  21 
B  22 
8  ?3 

8  24 

8  25- 
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SUBROUTINE  CDSF  (CD. AM. ANR) 

COMPUTES  INVISCIO  DRAG 

COMMON  /8/  AREA. RAD. SLANT. SNR. THETC.V 

COMMON  /C/  CPM.CP.CAN.CP9.C0F 

PI*3. 14159 

G«32.174 

AN0SE»PI«(SNR*«2j 
WAREA»PI*RAO»SLANT 
AHACHsAM 
PR* .89 

CPM*1 .42857/ (AM**2) • { (6.* (AM**2> /5. > **3.S« (6./ (7.* (Am**?) -1 . ) ) ••?. 
1 5- 1 . ) 

CP*4.*(SIN(THETcn**2*(2.5*8.*AMACH*SIN(THETC))/(l.*16.*AMACH*SlN( 
1 THETC ) } 

CAN* ( ,55-CP) *ANOSE/AREA 
CPB«1 .001/AM**? 

CALL  COFM  (COF) 

CD*CP*C»N*CPR*C0F 

RETURN 

END 


C 

C 

C 

r- 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


1 

2 

3 

4 

5 

6 
7 

e 

9 

10 

n 

12 

13 

14 

15 

16 

17 

18 
19 
20- 
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SUBROUTINE  LONG  (DLG.XE.YE.PI) 

c  THIS  SUB  COMPUTES  LONGITUDE  0  1 

IE  (YE)  1*6*4  0  2 

1  IE  (XE)  2.9*3  D  3 

?  DL6»-PI*ATXN(YE/XE)  0  ♦ 

RETURN  0  5 

3  OLG.XTXNCYE/xE)  0  6 

RETURN  0  7 

IE  (XE)  5.10.3  D  8 

DLG»PI**TAN ( YE/xE)  0  9 

RETURN  D  13 

IE  (XE)  7.7.8  D  11 

0L6a.PI  0  12 

RETURN  D  13 

*  OLG-0.0  0  1* 

RETURN  0  15 

9  0L6— PI/?.  0  16 

RETURN  0  17 

10  OLG-PI/2.  D  18 

RETURN  D  19 

End  o  20 

D  21- 
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UU  V  u 


1 


2 


3 


SUBROUTINE  INITIAL  «SN»«TM£TCt5L*NT)  £ 

£ 

THIS  SUBROUTINE  PRCYlDES  f*£  NECESSARY  GEOMETRY  CALCULATIONS  *NO  E 
CHANGES  IN  Uw!T<  FOR  PROGRAM  STRAR6  E 

E 

COMMON  /D/  TN{SO.12j.I»0{5Otl2j#JIKT){12}.VCP(5Otl2>.vTC{56.12i.V»H  £ 

10CS0.12)  .0EL{50.12J  «0£ I50»  121  *S0QT  I 121  .XI 112)  .RS t!2)  .«L U2:  .SOLU2  £ 
2) »MT (125 tTHC12) .T«T(12J .SACP{50> .SATClSO) ♦SARH0153) •RNciai  £ 

COMMON  /F/  N0FX.AL6TH«R*0IfP0SCll)*A2N.6L*T  £ 

AZK.aZN/57.296  E 

6L*T*SLAT/5T.296  E 

AL6TH«P0SInC?x1/12.  E 

Im**C?X-1  E 

STH»SIN|Th£Tcj  E 

CTH«COS«Th£TC}  £ 

DO  1  1*1 • IM  £ 

RNCU  *6.6£«0T  E 

soLm*iPosu*n-Rosnn/ii2.G*cTH}  e 

as  c  1 3  *S«31*CTh»  I OOS 1 1 3  /12  .  O-SNR*  C 1 . 0-STH3  i  ‘STM/CT:*  E 

RNtNOFX3*5.0E*C6  E 

RN(NOFX«13*1.SE*06  E 

R5tN©FXl*SNR  E 

RSINOFX*1J*SnR  e 

sol  <nofx3  » ipos  m /12. e-s«a*  c 1 *c-sthi j /cth  e 

otOl*SN3*CTH« (AL8TH.SNR*I1.0«STM>  3«sTH/CTH  E 

I«im«0EX.2  E 

GO  2  1*3 »IM  E 

RLII3*RSII*13  £ 

RtlNOFX-13»RADI  £ 

RL!NQrx3«aSU1  £ 

GL*NT»0.0  E 

DO  3  I*l.MOFx  £ 

SLANT*SLAN7*S0Lf!)  £ 

RETURN  E 

EnO  E 


3 

2 

3 

* 

5 

6 

7 

8 
9 

10 

12 

12 

13 

1* 

15 

16 

17 

18 

19 

20 
21 
22 
23 
2* 

25 

26 

27 

28 
29 
39 

31 

32 

33 
3* 
15- 
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SUBROUTINE  CDFM  (CDF) 

J^UBDC0MPUT*S  skin  fiction  DRAG 
COMMON  /B/  AREA , RAD* SLANT  *  SNR » THETC, V 

^COM  .ON  /E/  nofxm,nal,alt,amach,qd.amache,tempe,rhoe,ve,visg,pbe,an 

TW/*»0  . 

CDF«0  »  0 

SR«SNR/RAO 

PI»3. 14159 

CALL  RHOF  (RHO.aLT) 

call  visco  (vis.ald 

CALL  VSO  (VS, ALT) 

CALL  PRESS  (PHE.ALT) 

CALL  TEMPA  (TEMP, ALT) 

AMACH«V/VS 

f  !si:4t::cr s  ,N  *  thetc  ”  **s  > 

»M»CHE.»«»CH.,VE«> • (TEMP/TEMPE)...i  . 

PI “1 716, 

am»amach 

am-amach 

u:«;?^r2''!,N,THETc,,"j--j . 

RHOE«PBE/ (R1*TEMPE) 

CALL  VI SCOT  (VISG,TEMr_> 

AREA»PI*RA0*«2 

ANRE»RHOE*VE«SLANT«32. 174/VIS6 
COTT-1 ./TAN(THETC) 

00  2  K«l,NOFXM 
DO  1  W*1.NAL 

IF  (DE(M,K) .LE.O.)  GO  TO  1 
TWA»TWA*TW(M,K>/8. 

CONTINUE 

CONTINUE 

RE«RHO*V«SLANT/VIS 

IF  (ANRE~2,0E*07)  3,4,4 

HSH«,5*,5*(TWA/TEMPE5*,0374*(AMACHE)**2 

i»i:SSH!^n,-88‘2*s"-“-M3S*s"**j‘7i>-383*sR««i-isa . . 

GO  TO  5 

HSH«,S*.5*(TwA/TEMPE)*,0388*(AMACHE)*#2 

U/HSH|,:"s9./!o??,3I'1VE/V’"1,8*IPBE'',,BE,"-9*,TE"F/TEH',E’"-5*-'1 
CDF«COFi*(l  (  .  8*  ,  052 -*AmACH)  »SR)  ) 

COF«COF 

return 

END 


F  10 
F  11 
F  12 
F  13 
F  14 
F  IS 
F  16 
F  17 
F  18 
F  19 
F  20 
F  21 
F  22 
F  23 
E  24 
F  25 
F  26 
F  27 
F  28 
F  29 
7  30 

F  31 
F  32 
F  33 
F  34 
F  35 
F  36 
F  37 
F  38 
F  39 
F  40 
F  41 
F  42 
F  43 
F  44 
F  45 
F  46 
F  47 
F  48 
F  49 
F  50 
F  51 
F  52- 


110 


SUBROUTINE  VI5COT  (VISG*T£Mp) 

THIS  SU8  COMPUTES  VISCOSITY  OP  AIR  VS  TEMP 

UNITS  ARE  LB/FT-ScC 

IP  ( TEMP-2460* )  2,1,1 

VISG»41 .0E-10»TEMP*2,4414E-0S 

RETURN 

IP  (TEMP-1160.)  4,3,3 

VISQ«100.?69E-10#T£MP*.97l08»'-05 

RETURN 

IP  (TEMP-760.)  fc,5»5 

VIS6«132.5£-10*T€MP<.,603£-05 

RETURN 

VISG«166.666E-1o#TEMP*.34334E-05 

return 

END 


6  10 
G  11 
6  12 
G  13 
G  14 
G  15- 


t 


SUBROUTINE  PRESS  (PE. ALT) 

C  COMPUTES  PRESSURE  PER  ALTITUDE 

C  LBS/ET2 

PE«2116.2«EXP(-*.25509E-05*ALT) 

RETURN 

END 
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SUBROUTINE  VSO  (VS, ALT) 

C  THIS  SUR  COMPUTES  SPEED  OF  SOUND 

C  UNITS  ARE  FT/SEC 

IP  (ALT-36000.)  1,1.2 

1  VS«1116.9-.M2E-02*ALT 
RETURN 

2  IP  (ALT-82000.)  3,2 ,4 

3  V<>968.5 
RETURN 

*  IP  (ALT-IS4000.)  5,5,6 

5  VS“81*.*«19lE-OP#ALT 
return 

6  IP  (ALT-172000.)  7,7,3 

7  VS-1106. 

RETURN 

o  IF  (ALT-262000.)  9,9,10 

9  VS»1527.-.245E-02»ALT 
RETURN 

10  VS«RA6. 

RETURN 

END 


I 

I 

1 

I 

I 

I 

I 

I 

I 

I 

I 

I 

I 

I 

I 

I 

I 

I 

1 

I 

I 


I 

| 
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SUBROUTINE  GCPB  (GCP.PQU 

C  THIS  SUB  COMPUTES  ABLATIVE  GAS  SPECIFIC  HEAT 
C  UNITS  APE  8TU/IB-DEG  R 

IF  1P8L-10.)  2*1.1 

1  GCP»- *228*AL0G (P8L) *3.22 
RETURN 

2  IF  {P8L-.009)  4,3.3 

3  GCP«-.30882*AL0G(PBU*3.4 
RETURN 

IF  (P8L-.0001)  6.5,5 
GCP».73T4*ALO0{PBU  *8.37 
RETURN 

IF  (PBL-1.0E-05)  8.7,7 
GCP».12333*ALOG(P8U *2.709 

return 

8  GCPal.284 

RETURN 
ENO 


J  1 
J  2 
J  3 
J  4 
J  5 
J  6 
J  7 
J  8 
J  *» 

J  10 

j  n 

J  12 
J  13 
J  14 
J  15 
J  18 
J  17 
J  18- 
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OO  *-«  f\J  4  l/>  O  K  CD  <> 


SUBROUTINE  TEMP*  (TEMP, ALT)  X  1 

THIS  SUB  COMPUTES  TEMP  OF  AIR  VS  ALT  TO  300*000  K  2 

TEMP  IS  IN  DEGREES  RANKINE  k  3 

IF  (ALT-261000.)  2*1.1  *  4 

TEMP-298.2  *  5 

RETURN  X  8 

IF  (ALT-175500.)  4,3.3  K  1 

TEMP«-2.463£-03*ALY*941.04  X  8 

RETURN  X  9 

IF  (ALT-154000. )  6,5.5  X  10 

TEMP-508.79  X  11 

RETURN  X  12 

IF  (ALT-85000.)  8,7.7  X  13 

TEMP-1. 627E-03*ALT*256. 05  X  14 

RETURN  X  15 

IF  (ALT-35700.)  10.9,9  X  16 

TEMP-389.99  X  17 

RETURN  X  18 

10  TEMP— 36.77E-04-ALT-518.7  X  19 

RETURN  X  20 

ENO  X  21- 
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* 


SUBROUTINE  TC*T  (TCA.TEMP)  L  j 

this  sub  computes  thermal  cond.  of  air  vs.  temp  l 

units  are  btu/sec-ft-degree  rankine  l 

IF  (TEMP-2960.)  2.1.1  L 

TCA»6.0E-06*TEMP*,0332*  L 

GO  TO  9  L 

IF  (TEMP-2460.)  4.3.3  L 

TCA«7.8E-06*TEMP*.028  L 

GO  TO  9  L 

IF  (TEMP-1960.)  6.5.5  L 

TCA-14.2E-06*TEhP*. 01217  L 

GO  TO  9  L 

IF  (TEMP-1160,)  8.7,7  L 

TCA*16.5E-06*TEmP.. 00766  L 

GO  TO  9  • 

TCA«19.28E-06*TEMP*.00AA4  L 

TCA«TCA/3600.  L 

RETijRN  l 

ENO  L 
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SUBROUTINE  VCPPC  (VCP.TEMP)  M  1 

THIS  SUS  COMPUTES  VIRGIN  SPECIFIC  HEAT  VS.  TEMP  M  2 

UNITS  ARE  BTu/LR-OEG  R  M3 

IF  ( TEMP-400G . J  2.1.1  M  * 

VCP».55E-04*TEMP*1.23  N  5 

RETURN  M  6 

IF  (TEMP-2000.)  4.3.3  M  7 

VCP«.!9E-03*TEMP*.69  M  8 

RETURN  M  9 

IF  (TEMP-1060.)  6.S.5  M  10 

VCP«. 554E-03*TEmP-. 0423  M  11 

RETURN  M  12 

IF  (TEMP-860.)  8.7.7  M  13 

VCP«.975E-03»TEmP-.486S  m  14 

RETURN  M  15 

IF  (TEMP-660.)  10-9,9  M  16 

VCP«.21£-03*TEMP*.1714  M  17 

RETuPN  M  18 

in  VCP«.375E-03»TEmP*.062S  m  19 

RETURN  m  20 

END  M  21- 
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4  CD  *4  <*  in  * 


subroutine  rhof  (Rho»alt>  n  i 

C  This  SUB  COMPUTES  OENSITy  of  *ir  vs.  altitude  n 

COMMON  /A/  RE  N 

6ALT*RE*ALT/ CRE*ALT>  N 

RMOOa. 0023769  N 

Ra53.352  N 

E«2.7183  N 

IF  (GALT-36089.)  1,1.2  N 

1  T»S18. 688-3. 56616*GALT/1000.  N 

RHOaRHOO* (T/518.668) **4,256  N 

RETURN  N 

2  IF  (GALT-82021.)  3,3.4  N 

3  EP2* (360»9.-GALT) / (R*389.98S)  N 

RHOaRHOO* .297069*E**EP2  N 

RETURN  N 

IF  (GALT-154199.)  5.5.6  * 

Ts254.98BM.A4S92<*GALT/lOOO.  N 

RHOaRHOO*. 0326657* ( T/389 . 988 ) *♦- 12.388  N 

RETURN  n 

IF  (GALT-173885.)  7,7.8  N 

EP4a (154199. -GALT)/(R*508. 788)  N 

RHOaRHOO*. 00l21l79*E**£P4  N 

RETURN  N 

IF  (GALT-259186.)  9.9,10  N 

T*938. 088-2. 46888*GALT/1000.  N 

RHOaRMOO*5.86784E-04a (T/508.788) **6.592  N 

RETURN  N 

10  If  (GALT-295276.)  11,11.12  N 

11  EP6a(259186.-GALT)/(H*298.186)  N 

RHOaRHOO*i.73274E-05*E**EP6  N 

RETURN  N 

12  Ta-349, 812*2. 19*56*GALT/1000.  N 

RHOaRHOO* 1 . 792683E-06* ( T/?9S .188) **-9 .5*1  N 

RETURN  N 

END  N 
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o  00  -4  o.  gn  ►  w  i\>  •—  on 


SUBROUTINE  VISCO  (YIS.ALT)  0  1 

this  sub  computes  viscosity  of  air  vs.  altitude  o  2 

UNITS  ARE  in  LB-SEC/FT**2  0  3 

IF  {ALT-270000.1  2.1.1  0  * 

VIS»2.3SE-07  0  5 

RETURN  0  6 

IF  {ALT-175500. l  *.3.3  0  7 

VIS«-1.509E-12*ALT*6.333E-07  0  8 

RETURN  0  9 

IF  {ALT-156000. 1  6.5.5  O  10 

VfS«3.68E-07  0  11 

RETURN  0  12 

IF  {ALT-83250. )  8.7,7  0  13 

VIS*.973E-12*ALT-2,172E-07  0  1* 

RETURN  O  is 

IF  (ALT-34500.)  10.9,9  0  16 

VIS-2.969E-07  0  17 

RETURN  0  18 

VIS»-2.12*E-12*ALTO.737e-07  0  19 

RETURN  0  20 

END  0  21- 
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t 


SUBROUTINE  CCPP  CCCP.TEND)  p  j 

This  SUR  CO^puTfS  CHAR  SPECIFIC  HEAT  vs.  temp  p  2 

UNITS  are  STu/LR-DEG  R  P  3 

IF  tT£MP-6000.)  2.1.1  p  * 

CCP«.4  p  5 

RETURN  p  5 

IF  (TENP-4000. 1  4.3.3  p  7 

CCP*.12St-04*TE«P*,3215  p  g 

RETiiRn  p  9 

IF  { TEMP-3000.)  6.5,5  p  io 

CCP*.30r£-0*«TE*<P*.2495  p  jj 

RETURN  p  ip 

IF  {TEMP-2000. )  8.7,7  p  13 

CCP*.47£— 04*T£MP»« ’  *  p  j* 

RETURN  p  i5 

IF  {TEMP-1060. )  10.9,9  p  i6 

CCP«.207£-03*T£wP-.l*  p 

RETURN  p  ip 

10  CCP*.96E-04*TEMo-. 02236  p  19 

RETURN  p  20 

END  p  21. 
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WMm 


SUBROUTINE  SSLP  !VRZ4.VRX4*B£T*PI1 
C  COMPUTES  ANGLE  OF  SIDESLIP 

IF  {VSZ*>  9*5*1 

1  IF  (VRX4)  4*3*2 

2  3ET-ATANCVRX4/VRZ4) 

RETURN 

3  9ET-0. 

RETURN 

4  3ET-ATAN{ VRX4/VRZ4) 

RETURN 

6  IF  (VRX41  6*7*6 

*■  8ET--PI/2* 

RETURN 

7  BET-0. 

RETURN 

8  BET-Pl/2. 
return 

9  IF  (VRX4)  10.11,12 

1*5  3ET»PI«*T4N{VRX4/VRZ4} 

RETURN 
11  8ET.PI 

return 

1?  BET-P1*ATAN(VRX4/VRZ41 

return 

END 


0  1 
a  2 

0  3 

O  4 
0  S 
Q  6 
Q  7 

o  e 
a  9 
a  20 
o  n 
a  12 
a  13 

a  14 

a  is 
0  16 
a  it 

C  IB 

a  19 
a  20 
a  2i 
a  22 
a  23 

a  24 
a  25- 


SUBROUTINE  VTCPC  (VTC.TEMP) 
UNITS  ARE  BTU/FT-SEC-OEO  R 
IP  (TEMP-4000.)  2*1*1 
VTC*.705E-07*TEmP*,58£-04 

return 

IF  (TEMP-2000.)  4.3»3 

VTC«.69S£-07*TEmP*.62E-04 

RETURN 

VTC*.7376E-07*TE«P*.5347E-04 

return 

END 


R  8 
R  9 
R  10 
R  11* 
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SUBROUTINE  ACMQ  (C,A) 

C  COMPUTES  DAMPING  COEFFICIENT 

IF  (A. ST. 17.)  SO  TO  1 
IF  (A.GT.9.)  SO  TO  2 

IF  (A.GT.T.)  SO  TO  3 

IF  (A.GT.S.)  GO  TO  4 

IF  (A.GT .4. )  GO  TO  5 

IF  (A.GT. 3.)  GO  TO  6 

C«2.5*A-10. 

RETURN 

1  C--2.4 

return 

2  C«.fl5*A«3.2S 
RETURN 

3  C«.l*A-3.7 

return 

*  C*.3*A«5.3 

RETURN 

5  C«.7*A-?.3 
RETURN 

6  C«A-8.5 
return 
END 


S  1 
S  2 
S  3 
S  4 
S  5 
S  6 
S  7 
S  8 
$  9 

S  10 
S  11 
S  12 
S  13 
S  14 
S  15 
S  16 
S  17 
S  18 
S  19 
S  20 
S  21 
S  22 
S  23- 


SUBROUTINE  GASTC  <GTC»GCP} 

C  COMPUTES  ABLATIVE  GAS  THERMAL 

C  CONDUCTIVITY.  8TU/SEC-FT«D£G.R. 

GTC«A,*29E-05« (GCP*«?97) 

RETURN 

END 
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o  o 


SUBROUTINE  ACNA  (C,AJ  U  1 

COMPUTES  NORMAL  FORrE  COEFF.  U  2 

CNA  PER  DEGREE  U  3 

C«.03*  U  4 

RETURN  U  5 

END  U  6- 
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SUBROUTINE  CTCPC  (CTC.TENP) 

c  this  sub  computes  char  Thermal  conductivity  ve*  temp 

c  units  are  btu/ft-sec-oeq  R 

IP  < TEMP-6000*)  1*2,2 

CTC»,56£-07*T£MP*.053E-03 

return 

2  CTC*.389E-03 

return 

END 


V  1 

V  2 

V  3 

V  A 

V  5 

V  6 

V  7 

V  8 

V  9- 
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o  o 


SUBROUTINE  CPATJ  (S.T) 

COMPUTES  SPEC  MEAT  Of  ATj  GRAPHITE  VS.  TE«P 
UNITS  BTU/L8-OE6-R 
IF  (T-2000.)  1,2.2 

1  IF  JT-1000.)  4,3,3 

2  S»3,482E-05*T«. 35536 
RETURN 

3  S«1.25E-04«T*.l75 
return 

4  S»3.4O9E-O4»T-.O409 
return 

END 


127 


-*> 


I 


SUBROUTINE  TCATJ  IC.T) 

C  CONFUTES  THER.  CONO.  OF  ATJ  GRAPHITE  VS.  TEMP 
IF  IT-3460.)  1*3.3 

1  IF  IT-2460.)  2.4.4 

2  IF  IT-1460.)  6.5.5 

3  C*t2.E-0S»T*.28o8>/60. 

RETURN 

4  Ca.35/60. 

RETURN 

5  C* (-1 «84E-0**T*. 80264) /50, 

RETURN 

6  C«C-3.82E-04*T*1.091)/60. 

return 

END 


X  I 
X  2 
X  3 
X  4 
X  b 
X  6 
X  7 
X  8 
X  9 
X  10 
X  11 
X  12 
X  13 
X  14- 
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ATItvoo  (In tui  cnrilH  mm  bo  rumovtv]  mid  rtiml  Into  proftrnm  by  multl-cmm  rood  ritntvnmiitn. 


PRINT  87,  (SD0T(I),I=1,N0FXP)  A  773 
PRINT  70,  ((TW(I, J) , J>=l,NOFX) ,I=1,NP2)  A  774 
PRINT  88,  TW(NALCl.NOFX) ,TW(NALC2,NOFX)  A  775 
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FORMAT  (1H1)  A  784 
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APPENDIX  III 
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